PRÁCTICA ESTADÍSTICA DESCRIPTIVA, REGRESIÓN LINEAL MÚLTIPLE Y REMUESTREO.
EJERCICIO 1.
1.A.
#install.packages("ISLR")
library(ISLR)
data(Credit)summary(Credit)## ID Income Limit Rating
## Min. : 1.0 Min. : 10.35 Min. : 855 Min. : 93.0
## 1st Qu.:100.8 1st Qu.: 21.01 1st Qu.: 3088 1st Qu.:247.2
## Median :200.5 Median : 33.12 Median : 4622 Median :344.0
## Mean :200.5 Mean : 45.22 Mean : 4736 Mean :354.9
## 3rd Qu.:300.2 3rd Qu.: 57.47 3rd Qu.: 5873 3rd Qu.:437.2
## Max. :400.0 Max. :186.63 Max. :13913 Max. :982.0
## Cards Age Education Gender Student
## Min. :1.000 Min. :23.00 Min. : 5.00 Male :193 No :360
## 1st Qu.:2.000 1st Qu.:41.75 1st Qu.:11.00 Female:207 Yes: 40
## Median :3.000 Median :56.00 Median :14.00
## Mean :2.958 Mean :55.67 Mean :13.45
## 3rd Qu.:4.000 3rd Qu.:70.00 3rd Qu.:16.00
## Max. :9.000 Max. :98.00 Max. :20.00
## Married Ethnicity Balance
## No :155 African American: 99 Min. : 0.00
## Yes:245 Asian :102 1st Qu.: 68.75
## Caucasian :199 Median : 459.50
## Mean : 520.01
## 3rd Qu.: 863.00
## Max. :1999.00En este primer resumen de los datos podemos ver que tenemos algunas variables continuas, como Income, Limit, Rating (Credit Rating) y Balance. Por otro lado, tenemos variables discretas como el ID (No útil, sólo sirve para identificar registros), Cards (Número de tarjetas de crédito), Age (al estar puesta en número de años, y no poder tomar cualquier valor, sino solo valores enteros, es discreta), Education (años de educación). Además, tenemos variables categóricas, como es el caso de las binarias Gender, Student y Married; por otro lado tenemos la categórica de tres niveles Ethnicity.
Para utilizar el qqPlot que nos dé además un intervalo de confianza para la desviación de las colas, y poder así entender un poco más de nuestra distribución en cada una de las variables del set de datos, utilizaremos la librería car.
También usaremos la librería fitdistrplus, que nos ayudará a entender la distribución de nuestros datos.
#install.packages("car")
#install.packages("fitdistrplus")
if(!require(fitdistrplus)) {
install.packages("fitdistrplus")
}## Loading required package: fitdistrplus## Loading required package: MASS## Loading required package: survival## Loading required package: npsurv## Loading required package: lseiif(!require(car)) {
install.packages("car")
}## Loading required package: car## Loading required package: carDatalibrary(car)
library(fitdistrplus)credit <- data.frame(Credit)
credit$ID <- as.character(credit$ID)
#la transformamos en caracter; aunque esta variable no va a ser de utilidad.
boxplot(credit$Income, col = "blue")
hist(credit$Income, col = "grey")
qqPlot(credit$Income)
## [1] 29 324descdist(credit$Income, discrete = FALSE)
## summary statistics
## ------
## min: 10.354 max: 186.634
## median: 33.1155
## mean: 45.21889
## estimated sd: 35.24427
## estimated skewness: 1.742117
## estimated kurtosis: 5.947476#aquí, además, podemos ver los estadísticos descriptivos básicos de la variable, como la mediana, la media, y la estimación de la desviación estándar, la asimetría y la curtosis.
#ks.test(credit$Income, pnorm)
if(!require(nortest)) {
install.packages("nortest")
}## Loading required package: nortestlibrary(nortest)
#hacemos un test Lilliefors
lillie.test(credit$Income)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: credit$Income
## D = 0.16612, p-value < 2.2e-16#también hacemos un test de Anderson-Darling
ad.test(credit$Income)##
## Anderson-Darling normality test
##
## data: credit$Income
## A = 23.011, p-value < 2.2e-16Como podemos ver en los dos diagramas de arriba, el Income de las observaciones tiene una cantidad importante de outliers; en el boxplot se observa claramente una asimetría de la distribución, que confirmamos al dibujar el histograma. Desde luego el Income no se distribuye como una Normal.
Como podemos ver en el qqPlot() de la librería car, ni siquiera se asemeja en lo más mínimo a una normal. Esto es comprensible pues el Income tiene un tope en 0, sin embargo no existe un máximo teórico para esta variable, ya que puede tomar prácticamente cualquier valor positivo (la pregunta ¿cuál es el máximo ingreso que una persona puede tener? no se puede contestar). Además, como bien explican Brynjolfsson & McAffee(2014) en The Second Machine Age, entre otros, la riqueza, y por lo tanto también el Income, no se distribuyen de forma normal; pues un porcentaje pequeño de la muestra acapara un porcentaje grande de la riqueza.
Vemos por el output de la función descdist() que la distribución, de las incluidas en este paquete, que mejor se ajusta a la variable Income es la beta; esto se debe a una curtosis (concentración de los datos al rededor de la media) más alta que en la normal, y una square skewness muy superior. Los causantes de estas magnitudes son los outliers que tenemos en la muestra. Cabe mencionar que hablamos de outliers en lo que sería una distribución normal, en una distribución de otro tipo como la que probablemente tenga Income a juzgar por los datos es posible que estos valores no se consideren outliers.
No incluimos en el código el test de Kolgomorov-Smirnov, ya que tiene el problema de que existen valores duplicados en la muestra para la variable Income, y esto da un error. Por ello vamos a utilizar otros contrastes algo más robustos de cara a testar la normalidad de esta variable. En el Lilliefors podemos comprobar que efectivamente la variable Income está lejos de comportarse como una normal, ya que \(p-val < \alpha, \alpha = 0.05\). El Anderson-Darling test que hacemos a continuación también nos confirma esta última afirmación, por el mismo motivo.
boxplot(credit$Limit, col = "red")
library(ggplot2)
media <- mean(credit$Limit, na.rm = T)
std <- sd(credit$Limit, na.rm = T)
set.seed(1)
binwidth = 0.3
valores <- data.frame( valores = rnorm(n = nrow(credit), mean = media, sd = std))
credit$valores_aleatorios <- valores$valores
plt <- ggplot(data = credit, aes(x = Limit)) +
geom_histogram(bins = as.integer(sqrt(nrow(credit))), colour = "black",fill = "green" ,
aes(y = ..density.., fill = ..count..)) +
scale_fill_gradient("Count", low = "#DCDCDC", high="#7C7C7C") +
stat_function(fun = dnorm,
color = "blue",
args = list(mean = media,
sd = std)) +
geom_vline(xintercept = mean(na.omit(credit$Limit)), color = "red", show.legend = T) +
geom_vline(xintercept = median(na.omit(credit$Limit)), color = "blue", show.legend = T)
plt
qqPlot(credit$Limit) 
## [1] 324 29descdist(credit$Limit, discrete = F)
## summary statistics
## ------
## min: 855 max: 13913
## median: 4622.5
## mean: 4735.6
## estimated sd: 2308.199
## estimated skewness: 0.8374929
## estimated kurtosis: 4.004247lillie.test(credit$Limit)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: credit$Limit
## D = 0.072122, p-value = 3.454e-05ad.test(credit$Limit)##
## Anderson-Darling normality test
##
## data: credit$Limit
## A = 3.0579, p-value = 1.101e-07De nuevo observamos una asimetría a la derecha para la variable Limit, como podemos ver tanto en el boxplot como en el histograma. Los valores de las colas no se dispersan tanto de la normal como para la variable Income, pero siguen estando fuera del intervalo de confianza al 95%. En el último de estos 4 gráficos podemos ver que la distribución que más se ajusta a esta variable es una gamma, quedando también cerca de la curtosis y asimetría típica de una lognormal.
Estos resultados tienen bastante sentido, ya que esta variable representa el límite de crédito del cliente, que esperamos que se distribuya de forma parecida al Income; con una salvedad, y es que para determinados Incomes altos el Límite de crédito no aumentará tan rápido como para diferencias menores de Income de unos clientes a otros.
Viendo los resultados del Lilliefors test y el Anderson-Darling test, podemos confirmar que la variable Limit tampoco se distribuye como una normal.
boxplot(credit$Rating, col = "blue")
hist(credit$Rating, col = "green")
qqPlot(credit$Rating)
## [1] 324 29descdist(credit$Rating, discrete = F)
## summary statistics
## ------
## min: 93 max: 982
## median: 344
## mean: 354.94
## estimated sd: 154.7241
## estimated skewness: 0.8653935
## estimated kurtosis: 4.060947lillie.test(credit$Rating)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: credit$Rating
## D = 0.07039, p-value = 6.221e-05ad.test(credit$Rating)##
## Anderson-Darling normality test
##
## data: credit$Rating
## A = 3.2516, p-value = 3.718e-08El rating también sigue una distribución asimétrica a la derecha. Esto queda claro cuando hacemos el boxplot, que nos indica que un importante número de observaciones serían outliers en una Normal, y más claro aún cuando vemos el histograma. Vemos que la distribución que mejor se ajusta a estos datos es de nuevo una distribución Gamma, siendo la curtosis y la asimetría de esta variable también similares a los que tendría una distribución lognormal; como la distribución Weibull está cerca de estas dos, podríamos decir que habría que considerar que la variable Rating se puede distribuir también como uan Weibull.
Tanto el Lilliefors como el Anderson-Darling test nos indican que esta variable tampoco sigue la distribución normal.
boxplot(credit$Cards, col = "blue")
hist(credit$Cards, col = "green")
descdist(credit$Cards, discrete = T)
## summary statistics
## ------
## min: 1 max: 9
## median: 3
## mean: 2.9575
## estimated sd: 1.371275
## estimated skewness: 0.7919283
## estimated kurtosis: 3.943953Como el número de tarjetas de crédito siempre es positivo, tenemos un historama en forma de escalera, de izquierda a derecha; gráficamente, podríamos decir que se asemeja a una Power Law, aunque quizás el “codo” no está lo suficientemente marcado. No debemos olvidar que esta variable es discreta y con un rango bastante limitado como vimos en el summary(), es decir, es una variable que puede tomar pocos valores (no hay mucha gente con más de 1-2 tarjetas de crédito). No realizamos ningún test en esta variable porque, como vemos en los gráficos, toma tan sólo 8 valores y está claro que no se comporta de la forma “normal”; gráficamente esto queda patente. Más adelante decidiremos si es mejor tener esta variable como discreta numérica o como categórica.
boxplot(credit$Age, col = "blue")
hist(credit$Age, col = "green")
qqPlot(credit$Age)
## [1] 324 210descdist(credit$Age, discrete = T)
## summary statistics
## ------
## min: 23 max: 98
## median: 56
## mean: 55.6675
## estimated sd: 17.24981
## estimated skewness: 0.01149591
## estimated kurtosis: 1.93442Sin embargo la variable Age sí parece comportarse más o menos como una normal, con más observaciones en el centro de la distribución que en las “colas”; al hacer el boxplot vemos que en esta variable no se observan outliers a priori, y que existe quizás una pequeña asimetría a la izquierda, ya que la mediana es algo mayor a la media (la diferencia es muy pequeña).
Sin embargo, tal y como observamos en el resultado de descdist(), a lo que más se asemeja esta variable es a una normal, con algo menos de curtosis (esto lo vemos en el histograma, los datos no están tan centrados al rededor de la media como en una normal), también queda claro al ver el qqPlot, ya que en las colas tenemos más observaciones de las esperadas; las colas de la distribución de la variable no coinciden con las de una distribución Normal.
boxplot(credit$Education, col = "blue")
hist(credit$Education, col = "green")
qqPlot(credit$Education)
## [1] 31 225descdist(credit$Education, discrete = T)
## summary statistics
## ------
## min: 5 max: 20
## median: 14
## mean: 13.45
## estimated sd: 3.125207
## estimated skewness: -0.3292118
## estimated kurtosis: 2.421879Vemos que la variable Education, que representa el número de años de educación, toma valores entre 5 y 20. La distribución es asimétrica a la izquierda, como podemos ver en el boxplot y en el posterior histograma. Hay poca gente con más de 15 años de educación. Sin embargo, tal y como vemos en el último gráfico, podríamos decir que pese a ser una variable discreta, toma suficientes valores, y estos están distribuidos de tal forma que siguen un comportamiento “normal”.
ggplot(data = credit, aes(x = factor(Gender))) +
geom_bar(fill = "orange")
Vemos que la proporción de hombres y mujeres es prácticamente la misma; con algunas más mujeres que hombres. La representación de ambas categorías o clases es suficiente en la muestra, lo cual es positivo si utilizamos esta variable para construir modelos.
ggplot(data = credit, aes(x = factor(Student))) +
geom_bar(fill = "coral")
Vemos que tenemos muy pocos estudiantes en la muestra; esto significa que la mayoría de los clientes del banco (al menos basándonos en esta muestra de 400 clientes que nos han dado) no son estudiantes. Es razonable, ya que pocos estudiantes piden un crédito a no ser que sea para un préstamo de estudios. Una reflexión que puede surgir a partir del gráfico posterior es que el banco tiene un perfil de cliente más de gente a partir de una edad (como vimos en el gráfico de Age), y trabajando o jubilados. Esto nos da una idea del tipo de clientes que tiene el banco.
ggplot(data = credit, aes(x = factor(Married))) +
geom_bar(fill = "lightblue")
Además, gran parte de los clientes del banco (más de la mitad), están casados. Normalmente se asocia el hecho de estar casado a una mayor fiabilidad a la hora de repagar los préstamos etc, a un mayor compromiso, por lo tanto tiene sentido que el banco tenga una mayoría de clientes casados. Además, a juzgar por los datos de la edad y los estudiantes, la mayoría de los clientes del banco trabajan y tienen una edad superior a 30 años, por lo tanto están la mayoría en edad de estar casados.
ggplot(data = credit, aes(x = factor(Ethnicity))) +
geom_bar(fill = "darkblue")
Viendo la distribución de las etnias a las que pertenecen los clientes, es muy probable que los datos provengan de EEUU, uno de los países con más mezcla étnica del mundo; especialmente con muchos asiáticos y afro-americanos. Esto explica que, pese a haber una mayoría de caucásicos, también tienen un importante número de clientes de estas dos etnias; aproximadamente 1/2 caucásicos y 1/4 cada una de las otras dos etnias. Esto es positivo, ya que en la muestra tenemos suficientes casos de cada uno de los 3 niveles de este factor.
boxplot(credit$Balance, col = "blue")
media = mean(credit$Balance, na.rm = T)
std = sd(credit$Balance, na.rm = T)
set.seed(1)
binwidth = 0.3
#hist(credit$Balance, col = "green")
valores = data.frame( valores = rnorm(n = nrow(credit), mean = media, sd = std))
credit$valores_aleatorios <- valores$valores
plt <- ggplot(data = credit, aes(x = Balance)) +
geom_histogram(bins = as.integer(sqrt(nrow(credit))), colour = "black",fill = "green" ,
aes(y = ..density.., fill = ..count..)) +
scale_fill_gradient("Count", low = "#DCDCDC", high="#7C7C7C") +
stat_function(fun = dnorm,
color = "blue",
args = list(mean = media,
sd = std)) +
geom_vline(xintercept = mean(na.omit(credit$Balance)), color = "red", show.legend = T) +
geom_vline(xintercept = median(na.omit(credit$Balance)), color = "blue", show.legend = T)
plt
qqPlot(credit$Balance)
## [1] 324 29descdist(credit$Balance, discrete = F)
## summary statistics
## ------
## min: 0 max: 1999
## median: 459.5
## mean: 520.015
## estimated sd: 459.7589
## estimated skewness: 0.5845951
## estimated kurtosis: 2.472142lillie.test(credit$Balance)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: credit$Balance
## D = 0.12901, p-value < 2.2e-16ad.test(credit$Balance)##
## Anderson-Darling normality test
##
## data: credit$Balance
## A = 9.3916, p-value < 2.2e-16Vemos que la variable Balance en ningún caso se distribuye parecido a una normal. Tiene un mínimo en 0, en el que encontramos una gran cantidad de observaciones. Lo que vemos en esta variable es que todos los clientes que han cogido tienen mínimo un crédito de 0, es decir, que no hay ningún cliente que no tenga al menos tanto dinero en la cuenta de media que el que utiliza para los pagos (cuando el Balance es positivo, significa que el cliente ha tenido de media ese crédito en exceso ($) para el período).
Vemos que las variables Limit y Balance se distribuyen de forma relativamente parecida, tenemos varias variables distribuidas como beta, lognormal o gamma (o Weibull); esto sería de gran ayuda para que la distribución multidimensional al realizar el modelo de regresión posterior fuera una gaussiana en p dimensiones (dependiendo de las p variables que incluyamos en el modelo); de esta forma se aseguraría con mayor certeza el cumplimiento de los supuestos del General Linear Model, y el mejor estimador estimador lineal de las \(\beta\) sería el OLS.
Con el Lilliefors test y el Anderson-Darling test podemos comprobar, viendo el p-value, que la variable Balance no se comporta como una normal (de una forma más completa y correcta: dada la muestra de datos estudiada, y según los tests mencionados previamente, tenemos suficientes evidencias para rechazar la Hipótesis Nula de Normalidad en la variable Balance).
sum(is.na(credit))## [1] 0credit$valores_aleatorios <- NULL #esta variable la eliminamos ya, sólo la necesitábamos para hacer un gráfico. Vemos que no tenemos ningún valor ausente; por lo tanto podemos proceder a continuar con el ejercicio sin necesidad de imputar valores.
library(GGally)
ggcorr(cor(credit[ , unlist(lapply(credit, is.numeric))]), label = T)
cor(credit$Balance, credit$Limit)## [1] 0.8616973cor(credit$Balance, credit$Rating)## [1] 0.8636252En este gráfico de correlaciones vemos que la variable Balance sólo tiene una correlación significativa, en ambos casos positiva, con Limit (0.86) y Rating (0.86); también tiene correlación positiva, aunque algo más suave, con Income. Income, a su vez, está bastante relacionada con Limit y Rating. Debemos tener en cuenta las relaciones lineales observadas en esta matriz de correlación para la construcción del modelo, en consideraciones tales como la multicolinealidad.
1.B.
Comenzamos con el modelo simple inicial:
credit$Married <- as.factor(credit$Married)
#inicialmente vamos a incluir Cards como una variable categórica.
credit$Cards <- as.factor(credit$Cards)
linear_reg1 <- lm(Balance ~ Limit + Cards + Married, data = credit)
summary(linear_reg1)##
## Call:
## lm(formula = Balance ~ Limit + Cards + Married, data = credit)
##
## Residuals:
## Min 1Q Median 3Q Max
## -650.10 -132.62 -9.72 132.83 746.11
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.945e+02 4.414e+01 -6.673 8.66e-11 ***
## Limit 1.725e-01 5.093e-03 33.868 < 2e-16 ***
## Cards2 -1.918e+01 3.907e+01 -0.491 0.624
## Cards3 -1.026e+00 3.936e+01 -0.026 0.979
## Cards4 5.189e+01 4.243e+01 1.223 0.222
## Cards5 1.165e+02 5.142e+01 2.266 0.024 *
## Cards6 9.019e+01 7.725e+01 1.168 0.244
## Cards7 8.151e+00 1.213e+02 0.067 0.946
## Cards8 2.636e+02 2.345e+02 1.124 0.262
## Cards9 1.434e+02 2.341e+02 0.613 0.541
## MarriedYes -3.166e+01 2.407e+01 -1.315 0.189
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 231.6 on 389 degrees of freedom
## Multiple R-squared: 0.7526, Adjusted R-squared: 0.7462
## F-statistic: 118.3 on 10 and 389 DF, p-value: < 2.2e-161.C.
Modelo teórico:
\[ Balance_i = \hat{\beta_0} + \hat{\beta_1}*Limit_i + \hat{\beta_2}*Cards_2i + \hat{\beta_3}*Cards_3i + \hat{\beta_4}*Cards_4i + \hat{\beta_5}*Cards_5i + \hat{\beta_6}*Cards_6i + \hat{\beta_7}*Cards_7i + \hat{\beta_8}*Cards_8i + \hat{\beta_9}*Cards_9i + \hat{\beta_10}*MarriedYes_i + \hat{\epsilon_i} \]
Modelo numérico:
\[ Balance_i = -294.5 + 0.1725*Limit_i + -19.18*Cards_2i + -1.026*Cards_3i + 51.89*Cards_4i + 116.5*Cards_5i + 90.19*Cards_6i + 8.151*Cards_7i + 263.6*Cards_8i + 143.4*Cards_9i + -31.66*MarriedYes_i + \hat{\epsilon_i} \] Interpretación de Coeficientes:
En primer lugar, dados los p-value de las variables independientes utilizadas en el modelo, podemos decir que tan sólo las variables Intercept (La constante), Limit y Cards5 (para la cual tendremos que crear una nueva variable e incluir únicamente esta en el modelo, en lugar de todos los niveles de tarjetas de crédito) son significativas. El resto de variables introducidas nos salen no significativas.
Pese a que antes de poder analizar profundamente los coeficientes significativos del modelo hay que eliminar las variables no signficativas y quedarse sólo con las que tengan un p-valor inferior a 0.05 (\(\alpha = 0.5\)) (una vez hecho esto los coeficientes del modelo cambiarán), procedemos a interpretar los coeficientes actuales:
–> Limit : Manteniendo todo el resto de variables constantes, un incremento en 1 unidad de limit está relacionado con un incremento de 0.1725 de Balance. Esto tiene sentido, ya que aquellos clientes que mantienen un balance mayor en la tarjeta de crédito tendrán un límite mayor de crédito, son más fiables.
–> Cards5: Manteniendo todo el resto de variables constantes, el hecho de que el individuo i tenga 5 tarjetas de crédito está asociado con un Balance 116.5 unidades mayor con respecto a si tuviera 1, 2, 3, 4, 6. 7 u 8 tarjetas de crédito.
El resto de coeficientes no tendría sentido interpretarlos debido a que no son estadísticamente significativos en la muestra de 400 observaciones con la que estamos trabajando. Para las variables Cards_NumberOfCards la interpretación seguiría la misma lógica que para la variable Cards5 incluida arriba. De igual forma, la variable MarriedYes se interpretaría como que, manteniendo todo el resto de variables constantes, estar casado estaría asociado con un Balance -31.66 unidades menor, de media, con respecto a no estarlo. Sin embargo, al no ser estadísticamente significativo (es 0), debemos decir que estar casado no tiene ningún efecto sobre el Balance de los clientes.
El \(R^2\) de este primer modelo es de 0.7526, mientras que el ajustado es de 0.742; esta diferencia no es excesiva, pero sí lo suficientemente como para poder afirmar lo mismo que nos llevan a afirmar los p-values, y es que al modelo hay que eliminarle variables, lo cual seguramente reduzca el \(R^2\), pero será un modelo más robusto y sólo compuesto de variables significativas.
hist(scale(credit$Balance), col = "blue")
minmaxscaler <- function(x){
(x-min(x))/(max(x)-min(x))
}
hist(minmaxscaler(credit$Balance), col = "lightblue")
hist(scale(credit$Limit), col = "green")
hist(minmaxscaler(credit$Limit), col = "lightgreen")
A juzgar por los gráficos, tendría más sentido, en caso de querer centrar nuestros datos y comprobar si esto tiene algún efecto positivo en la estimación del modelo, escalarlos con la función scale(), ya que gráficamente esta transformación hace que las variables Balance y Limit se parezcan algo más entre sí y a una Normal que la transformación con la función minmaxscaler()).
new_df <- credit[ , c("Balance", "Limit", "Married", "Cards")]
cards5num <- ifelse(as.numeric(new_df$Cards) == 5, 1, 0) #lo vamos a necesitar luego para hacer plots y demás, porque cuando le pasemos a esto de abajo as.numeric(Cards5) te devuelve 1 y 2, por los niveles que toma, en lugar de 0 y 1.
new_df$Cards5 <- as.factor(ifelse(as.numeric(new_df$Cards) == 5, 1, 0))
new_df$Cards <- NULL
linear_reg <- lm(Balance ~ Limit + Cards5 + Married, data = new_df)
summary(linear_reg)##
## Call:
## lm(formula = Balance ~ Limit + Cards5 + Married, data = new_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -659.46 -134.15 -12.59 137.86 763.57
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.872e+02 3.023e+01 -9.501 < 2e-16 ***
## Limit 1.728e-01 5.042e-03 34.278 < 2e-16 ***
## Cards51 1.090e+02 4.169e+01 2.615 0.00927 **
## MarriedYes -3.332e+01 2.381e+01 -1.399 0.16251
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 231.7 on 396 degrees of freedom
## Multiple R-squared: 0.7479, Adjusted R-squared: 0.746
## F-statistic: 391.7 on 3 and 396 DF, p-value: < 2.2e-161.D.
Vemos que en este modelo de nuevo MarriedYes sale no significativo, por lo tanto lo eliminamos.
linear_reg <- lm(Balance ~ Limit + Cards5, data = new_df)
summary(linear_reg)##
## Call:
## lm(formula = Balance ~ Limit + Cards5, data = new_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -671.62 -138.96 -9.27 139.19 784.12
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.063e+02 2.702e+01 -11.335 <2e-16 ***
## Limit 1.726e-01 5.045e-03 34.209 <2e-16 ***
## Cards51 1.065e+02 4.170e+01 2.555 0.011 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 232 on 397 degrees of freedom
## Multiple R-squared: 0.7467, Adjusted R-squared: 0.7454
## F-statistic: 585.1 on 2 and 397 DF, p-value: < 2.2e-16plot(x = 1:length(linear_reg$residuals), y = linear_reg$residuals, col = "red", main = "Residuals Plot", ylab="Residuals", xlab = "Observation")
abline(h = 3*sd(linear_reg$residuals), col = "blue")
abline(h = mean(linear_reg$residuals), col = "black")
abline(h = -3*sd(linear_reg$residuals), col = "blue")
En este modelo podemos ver que el \(R^2\) baja a 0.7467, pero el \(R^2adj\) es de 0.7454; podemos explicar, en teoría, el 74.54% de la varianza de la variable Balance con las variables Limit y Cards5, junto con la constante -306.3. Además, esta pequeña diferencia entre el \(R^2\) y el \(R^2adj\) significa que tenemos poca penalización por la inclusión de variables no significativas, por lo tanto nos reafianza la decisión de dejar como únicas variables del modelo Limit y Cards5.
Parte del código inferior para los gráficos en 3 dimensiones está basado en: http://www.sthda.com/english/wiki/a-complete-guide-to-3d-visualization-device-system-in-r-r-software-and-data-visualization
library(rgl)
library(knitr)
knit_hooks$set(webgl = hook_webgl)
#the below function was defined by http://www.sthda.com/english/wiki/a-complete-guide-to-3d-visualization-device-system-in-r-r-software-and-data-visualization
rgl_init <- function(new.device = FALSE, bg = "white", width = 640) {
if( new.device | rgl.cur() == 0 ) {
rgl.open()
par3d(windowRect = 50 + c( 0, 0, width, width ) )
rgl.bg(color = bg )
}
rgl.clear(type = c("shapes", "bboxdeco"))
rgl.viewpoint(theta = 15, phi = 20, zoom = 0.7)
}rgl_init()
rgl.spheres(x = scale(new_df$Balance), y = new_df$Cards5, z = scale(new_df$Limit), r = 0.1, color = "yellow") # Scatter plot
rgl.bbox(color = "#333377") # Add bounding box decorationEn este objeto se puede apreciar mucho mejor la relación tridimensional de las variables de estudio.
linear_reg <- lm(scale(Balance) ~ scale(Limit) + Cards5, data = new_df)
summary(linear_reg)##
## Call:
## lm(formula = scale(Balance) ~ scale(Limit) + Cards5, data = new_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.46081 -0.30225 -0.02016 0.30274 1.70550
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.01970 0.02638 -0.747 0.456
## scale(Limit) 0.86638 0.02533 34.209 <2e-16 ***
## Cards51 0.23172 0.09070 2.555 0.011 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5046 on 397 degrees of freedom
## Multiple R-squared: 0.7467, Adjusted R-squared: 0.7454
## F-statistic: 585.1 on 2 and 397 DF, p-value: < 2.2e-16library(plot3D)
scatter3D(z = scale(new_df$Limit), x = scale(new_df$Balance), y = as.numeric(new_df$Cards5), col = "blue")
scatter3D(z = scale(new_df$Limit), x = scale(new_df$Balance), y = as.numeric(new_df$Cards5), title = "3D Scatter Plot", add = F)
lines3D(x = linear_reg$coefficients[1]+linear_reg$coefficients[2]*scale(new_df$Limit) + linear_reg$coefficients[3]*as.numeric(new_df$Cards5), z = scale(new_df$Limit), y = as.numeric(new_df$Cards5), add = T)
Aprovechando que tenemos un modelo de pocas variables, he decidido graficarlo en 3 dimensiones; primero he hecho un scatter plot de las 3 variables, con los ejes de la siguiente forma:
Eje x: Balance escalado –> scale(Balance)
Eje y: Cards5
Eje z: Limit escalado –> scale(Limit)
Después he procedido a dibujar el mismo gráfico pero ajustando la recta de regresión que hemos sacado, de tal forma que podemos observar cómo quedan de dispersos nuestros datos al rededor de este plano (el que forma la recta que explica el eje x: Balance escalado, con los valores en cada uno de los ejes z:Limit escalado e y: Cards5). En caso de que tratemos Cards como una variable numérica, tendremos un plano, en caso contrario tenemos dos líneas de regresión.
Lo que ocurre es que al ser la variable Cards5 categórica binaria, tenemos dos líneas de regresión, la más oscura que es la de aquellos clientes que no tienen 5 tarjetas de crédito, y al rededor de la cual se aglutinan la mayor parte de los datos, y por otro lado la verde oscuro, ambas conectadas por líneas que las unen pues pertenecen ambas a x; esta es la recta resultante para aquellos clientes que sí tienen 5 tarjetas de crédito.
Al escalar las variables Balance y Limit el intercept pasa a ser 0. En este último modelo podemos ver que el \(R^2adj\) no se reduce.
1.E.
plot(linear_reg)
En estos plots podemos ver que ninguna observación supera la distancia de Cook, por lo que, pese a existir puntos con un leverage por encima de 0.03, siguiendo el criterio de la distancia de Cook no deberíamos quitar ninguna observación, por no considerarse que sesga el modelo lo suficiente.
Sin embargo, otra regla para eliminar observaciones utilizando la distancia de Cook, en lugar de cortar en 0.5 o en 0.1, es cortar cuando \(Cook's Distance > 4/n, n = number of observations\). Esto lo probaremos próximamente. Este método para eliminar observaciones utilizando la distancia de Cook viene muy bien explicado en este link.
library(car)
vif(linear_reg)## scale(Limit) Cards5
## 1.005269 1.005269#♦vemos que no hay colinealidad.
leveragePlots(linear_reg)
No tenemos problemas de colinealidad tampoco.
Ahora probaremos quitando los high leverage points de la forma descrita arriba:
HighLeveragePoints <- cooks.distance(linear_reg) > 4/nrow(new_df)
new_df_noleverage <- new_df[-HighLeveragePoints, ]
summary(lm(scale(Balance) ~ scale(Limit) + Cards5, data = new_df_noleverage))##
## Call:
## lm(formula = scale(Balance) ~ scale(Limit) + Cards5, data = new_df_noleverage)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.45926 -0.30519 -0.02019 0.30276 1.70379
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.01973 0.02642 -0.747 0.4556
## scale(Limit) 0.86635 0.02536 34.156 <2e-16 ***
## Cards51 0.23159 0.09073 2.552 0.0111 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5047 on 396 degrees of freedom
## Multiple R-squared: 0.7466, Adjusted R-squared: 0.7453
## F-statistic: 583.3 on 2 and 396 DF, p-value: < 2.2e-16Vemos que no mejora significativamente el ajuste de los datos tras eliminar los highleveragepoints.
par(mfrow=c(2,2)) # init 4 charts in 1 panel
plot(linear_reg)
Viendo los residuos frente a \(\hat{y}\) (\(\hat{Balance}\)), no observamos ninguna tendencia concreta. Esto es positivo ya que si observáramos un patrón común en toda esta nube de puntos, significaría que existe relación entre \(\hat{\epsilon}\) e \(\hat{y}\). Sin embargo, podemos ver que las “colas” de la distribución de los residuos se salen de lo que sería normal, como podemos ver en el Normal Q-Q plot (esquina superior derecha).
mean(linear_reg$residuals)## [1] 2.725196e-17cor(linear_reg$residuals, credit$Balance)## [1] 0.5033024cov(linear_reg$residuals, credit$Balance)## [1] 116.463fuera <- linear_reg$residuals[linear_reg$residuals > 2*sd(linear_reg$residuals) | linear_reg$residuals < -2*sd(linear_reg$residuals)]
n_fuera <- length(fuera)
print(n_fuera)## [1] 23plot(x = 1:length(linear_reg$residuals), y = linear_reg$residuals, col = "red")
abline(h = 2*sd(linear_reg$residuals), col = "blue")
abline(h = mean(linear_reg$residuals), col = "black")
abline(h = -2*sd(linear_reg$residuals), col = "blue")
Como podemos ver en los gráficos y en el número de valores de \(\hat{\epsilon}\) que se salen del intervalo al 99%, el efecto del escalado no es real, ya que tan sólo centra los datos al rededor de una escala más pequeña, pero se siguen quedando fuera estos valores (posiblemente atípicos) que tendremos que tratar posteriormente. El siguiente paso será realizar un test de homocedasticidad.
De todas formas, estamos omitiendo variables muy importantes para explicar el Balance seguro; esto provoca que se incumpla uno de los supuestos del Modelo Lineal General, y es que \[E[\epsilon] \neq 0\]. Esto se debe a que las variables omitidas se incluyen en el error; esto sesga el estimador OLS. Estimamos \(E[\hat{\epsilon}]\) y nos sale 0, por lo que la consideración teórica anterior podría no ser cierta. Sin embargo, como vemos en la correlación entre los residuos y la variable Balance, sí existe una cierta correlación entre estas dos variables; en teoría, \(E[(Balance - E[Balance])(\epsilon - E[\epsilon])] = 0\); sin embargo tanto la covarianza estimada (116.463) como la correlación estimada (0.50) nos indican que esta suposición es incorrecta, reforzando la afirmación de que estamos omitiendo importantes variables explicativas para el Balance.
if(!require(lmtest)) {
install.packages("lmtest")
}## Loading required package: lmtest## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericrequire(lmtest)
bptest(linear_reg)##
## studentized Breusch-Pagan test
##
## data: linear_reg
## BP = 18.378, df = 2, p-value = 0.0001021Como \(p-val < 0.05\), rechazamos la hipótesis nula de que existe homocedasticidad en el error del modelo.
Ahora vamos a probar quitando los valores que quedan fuera de \(\pm 3 \cdot \hat{\sigma_\epsilon}\), para ver qué tal funciona en ese caso el modelo.
take_outliers_out <- function(df, residuals = linear_reg$residuals) {
#Elimina aquellas variables que estén más allá de 3 sd (no hay ninguna observación por encima de 4sd)
if(is.data.frame(df)) {
new_df <- df[-which(residuals > 3*sd(residuals) | residuals < -3*sd(residuals)), ]
return(new_df)
} else {
new_vec <- df[-which(residuals > 3*sd(residuals) | residuals < -3*sd(residuals))]
return(new_vec)
}
}
new_df2 <- take_outliers_out(df = new_df)
model_complete_clean <- lm(Balance ~ Limit + Cards5, data = new_df2)
summary(model_complete_clean)##
## Call:
## lm(formula = Balance ~ Limit + Cards5, data = new_df2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -659.53 -134.05 -2.55 142.94 697.58
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.100e+02 2.568e+01 -12.071 < 2e-16 ***
## Limit 1.716e-01 4.795e-03 35.787 < 2e-16 ***
## Cards51 1.143e+02 3.963e+01 2.883 0.00415 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 220.4 on 393 degrees of freedom
## Multiple R-squared: 0.7652, Adjusted R-squared: 0.764
## F-statistic: 640.4 on 2 and 393 DF, p-value: < 2.2e-16shapiro.test(model_complete_clean$residuals)##
## Shapiro-Wilk normality test
##
## data: model_complete_clean$residuals
## W = 0.99141, p-value = 0.02144Como podemos ver, el modelo \(\hat{Balance_i} = - 310 + 0.1716*Limit_i + 114.3*Cards5_i\) tiene un \(R^2 = 0.7652\) y \(R^2adj = 0.764\); el ajuste mejora significativamente cuando quitamos esas observaciones.
Sin embargo, el Shapiro-Wilk test nos indica que los residuos no se comportan de forma normal, ya que \(0.02144 = p-value < \alpha = 0.05\). Por lo tanto se incumple uno de los supuestos del Teorema de Gauss-Markov; \(\epsilon \nsim N(0, \sigma^2_\epsilon)\). Esto implica que nuestro estimador OLS ya no será Best Linear Unbiased Estimator.
scatter3D(z = new_df$Limit, x = new_df$Balance, y = as.numeric(new_df$Cards5), col = "blue")
scatter3D(z = new_df2$Limit, x = new_df2$Balance, y = as.numeric(new_df2$Cards5), title = "3D Scatter Plot", add = F)
lines3D(x = model_complete_clean$coefficients[1]+model_complete_clean$coefficients[2]*new_df2$Limit + model_complete_clean$coefficients[3]*as.numeric(new_df2$Cards5), z = new_df2$Limit, y = as.numeric(new_df2$Cards5), add = T)
new_df2$Cards5 <- take_outliers_out(df = cards5num)
par(mfrow=c(2,1))
plot(x = new_df2$Limit[as.numeric(new_df2$Cards5) == 1], y = new_df2$Balance[as.numeric(new_df2$Cards5) == 1], col = "blue",
main = "Model for Cards5=1", xlab = "Limit", ylab = "Balance")
abline(-310 + 114.3, 0.1716)
plot(x = new_df2$Limit[as.numeric(new_df2$Cards5) == 0], y = new_df2$Balance[as.numeric(new_df2$Cards5) == 0], col = "red",
main = "Model for Cards5=0", xlab = "Limit", ylab = "Balance")
abline(-310, 0.1716)
Una de las conclusiones que se puede sacar del gráfico superior es que tenemos demasiadas pocas observaciones de clientes con cinco tarjetas de crédito en comparación con la muestra.
Un 8% es una proporción demasiado pequeña de una clase de una variable categórica, y por lo tanto no tenemos suficientes observaciones para estimar un modelo realmente robusto.
Por este motivo, vamos a proceder a realizar todo el análisis superior con un modelo en el que la variable Cards sea numérica en lugar de categórica. Además, probaremos a incluir pequeñas modificaciones en el modelo con el objetivo de mejorarlo. Un ejemplo de esto es incluir Limit^2.
credit$Cards <- as.numeric(credit$Cards)
#credit$Balance <- scale(credit$Balance)
#credit$Limit <- scale(credit$Limit)
modelo <- lm(Balance ~ Limit + I(Limit^2) + Cards, data = credit )
summary(modelo)##
## Call:
## lm(formula = Balance ~ Limit + I(Limit^2) + Cards, data = credit)
##
## Residuals:
## Min 1Q Median 3Q Max
## -664.09 -137.69 -9.32 136.52 773.48
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.882e+02 5.203e+01 -9.382 < 2e-16 ***
## Limit 2.226e-01 1.697e-02 13.122 < 2e-16 ***
## I(Limit^2) -4.490e-06 1.424e-06 -3.153 0.00174 **
## Cards 2.652e+01 8.346e+00 3.177 0.00160 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 228.6 on 396 degrees of freedom
## Multiple R-squared: 0.7547, Adjusted R-squared: 0.7528
## F-statistic: 406.1 on 3 and 396 DF, p-value: < 2.2e-16Como vemos, el modelo se ajusta algo mejor que cuando incluimos Cards5; es mejor incluir Cards como variable numérica.
plot(hatvalues(modelo), col = "blue")
En este gráfico vemos que hay unos 6-7 puntos con un hatvalue significativamente mayor al resto. Vamos a probar a eliminarlos con el fin de determinar si esto mejora la bondad de ajuste del modelo.
high_hatvals <- hatvalues(modelo)>0.05
credit_no_high_hatvals <- credit[!high_hatvals, ]
modelo_no_high_hatvals <- lm(Balance ~ Limit + I(Limit^2) + Cards, data = credit_no_high_hatvals)
summary(modelo_no_high_hatvals)##
## Call:
## lm(formula = Balance ~ Limit + I(Limit^2) + Cards, data = credit_no_high_hatvals)
##
## Residuals:
## Min 1Q Median 3Q Max
## -652.90 -139.90 -8.29 144.27 768.90
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.027e+02 5.784e+01 -8.692 < 2e-16 ***
## Limit 2.313e-01 2.128e-02 10.869 < 2e-16 ***
## I(Limit^2) -5.409e-06 1.992e-06 -2.716 0.00691 **
## Cards 2.593e+01 8.762e+00 2.959 0.00328 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 230.1 on 389 degrees of freedom
## Multiple R-squared: 0.7283, Adjusted R-squared: 0.7262
## F-statistic: 347.6 on 3 and 389 DF, p-value: < 2.2e-16plot(hatvalues(modelo_no_high_hatvals))
plot(modelo_no_high_hatvals)
La bondad de ajuste de hecho empeora cuando eliminamos esas observaciones, por lo que descartamos que quitar esas observaciones mejore el modelo; salen más hatvalues altos que en el anterior modelo, sugiriendo que los ajustes en el estimador tras el cambio en los datos sesgan al mismo, llevándole a predecir peor la variable Balance.
Vamos a quitarle los 3 valores que se salen de \(\pm 3*\sigma_\epsilon\), los cuales pueden ser considerados outliers.
df_no_outliers <- take_outliers_out(credit, residuals = modelo$residuals)
modelo2 <- lm(Balance ~ Limit + I(Limit^2) + Cards, data = df_no_outliers)
summary(modelo2)##
## Call:
## lm(formula = Balance ~ Limit + I(Limit^2) + Cards, data = df_no_outliers)
##
## Residuals:
## Min 1Q Median 3Q Max
## -656.03 -129.83 -4.99 141.58 682.23
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.839e+02 5.014e+01 -9.651 < 2e-16 ***
## Limit 2.156e-01 1.637e-02 13.170 < 2e-16 ***
## I(Limit^2) -3.964e-06 1.373e-06 -2.886 0.004115 **
## Cards 2.946e+01 8.079e+00 3.646 0.000302 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 220 on 393 degrees of freedom
## Multiple R-squared: 0.7677, Adjusted R-squared: 0.7659
## F-statistic: 432.8 on 3 and 393 DF, p-value: < 2.2e-16plot(modelo2)
Vemos que después de eliminar estas 3 observaciones mencionadas arriba el \(R^2\) aumenta, al igual que el ajustado. Podemos explicar el 76.59% de la varianza de Balance con este modelo:
\[ \hat{Balance} = -483.9 +0.2156 \cdot Limit - 3.964\cdot 10^-6 \cdot Limit^2 + 29.46 \cdot Cards \]
Vemos que un incremento de Limit en 1 unidad, manteniendo el resto de variables constantes, está asociado con un incremento de 0.2156 del Balance. El mismo tipo de razonamiento podría aplicarse a las otras dos variables; manteniéndose constantes el resto, el coeficiente estimado \(\beta\) de cada una de ellas tiene una relación de el tipo descrito arriba con la variable Balance.
Para finalizar, podríamos continuar con el análisis de este modelo comprobando, como para el modelo anterior, algo más exhaustivamente los leverage points, la homocedasticidad o la normalidad de los residuos.
vif(modelo2)## Limit I(Limit^2) Cards
## 11.760348 11.761499 1.000336shapiro.test(modelo2$residuals)##
## Shapiro-Wilk normality test
##
## data: modelo2$residuals
## W = 0.98707, p-value = 0.001335sum(cooks.distance(modelo2) > 4/nrow(df_no_outliers))## [1] 21bptest(modelo2)##
## studentized Breusch-Pagan test
##
## data: modelo2
## BP = 28.058, df = 3, p-value = 3.531e-06Podemos tener algo de colinealidad entre Limit y Limit^2, como era de esperar. Además, rechazamos la hipótesis de normalidad de los residuos y también la de homocedasticidad.
EJERCICIO 2
2.A.
library(mlbench)
data(BostonHousing)
housing <- data.frame(BostonHousing)
head(housing)## crim zn indus chas nox rm age dis rad tax ptratio b
## 1 0.00632 18 2.31 0 0.538 6.575 65.2 4.0900 1 296 15.3 396.90
## 2 0.02731 0 7.07 0 0.469 6.421 78.9 4.9671 2 242 17.8 396.90
## 3 0.02729 0 7.07 0 0.469 7.185 61.1 4.9671 2 242 17.8 392.83
## 4 0.03237 0 2.18 0 0.458 6.998 45.8 6.0622 3 222 18.7 394.63
## 5 0.06905 0 2.18 0 0.458 7.147 54.2 6.0622 3 222 18.7 396.90
## 6 0.02985 0 2.18 0 0.458 6.430 58.7 6.0622 3 222 18.7 394.12
## lstat medv
## 1 4.98 24.0
## 2 9.14 21.6
## 3 4.03 34.7
## 4 2.94 33.4
## 5 5.33 36.2
## 6 5.21 28.7summary(housing)## crim zn indus chas
## Min. : 0.00632 Min. : 0.00 Min. : 0.46 0:471
## 1st Qu.: 0.08204 1st Qu.: 0.00 1st Qu.: 5.19 1: 35
## Median : 0.25651 Median : 0.00 Median : 9.69
## Mean : 3.61352 Mean : 11.36 Mean :11.14
## 3rd Qu.: 3.67708 3rd Qu.: 12.50 3rd Qu.:18.10
## Max. :88.97620 Max. :100.00 Max. :27.74
## nox rm age dis
## Min. :0.3850 Min. :3.561 Min. : 2.90 Min. : 1.130
## 1st Qu.:0.4490 1st Qu.:5.886 1st Qu.: 45.02 1st Qu.: 2.100
## Median :0.5380 Median :6.208 Median : 77.50 Median : 3.207
## Mean :0.5547 Mean :6.285 Mean : 68.57 Mean : 3.795
## 3rd Qu.:0.6240 3rd Qu.:6.623 3rd Qu.: 94.08 3rd Qu.: 5.188
## Max. :0.8710 Max. :8.780 Max. :100.00 Max. :12.127
## rad tax ptratio b
## Min. : 1.000 Min. :187.0 Min. :12.60 Min. : 0.32
## 1st Qu.: 4.000 1st Qu.:279.0 1st Qu.:17.40 1st Qu.:375.38
## Median : 5.000 Median :330.0 Median :19.05 Median :391.44
## Mean : 9.549 Mean :408.2 Mean :18.46 Mean :356.67
## 3rd Qu.:24.000 3rd Qu.:666.0 3rd Qu.:20.20 3rd Qu.:396.23
## Max. :24.000 Max. :711.0 Max. :22.00 Max. :396.90
## lstat medv
## Min. : 1.73 Min. : 5.00
## 1st Qu.: 6.95 1st Qu.:17.02
## Median :11.36 Median :21.20
## Mean :12.65 Mean :22.53
## 3rd Qu.:16.95 3rd Qu.:25.00
## Max. :37.97 Max. :50.00Explicación Variables:
crim: per capita crime rate by town
zn: proportion of residential land zoned for lots over 25,000 sq.ft
indus: proportion of non-retail business acres per town
chas: Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)
nox: nitric oxides concentration (parts per 10 million)
rm: average number of rooms per dwelling
age: proportion of owner-occupied units built prior to 1940
dis: weighted distances to five Boston employment centres
rad: index of accessibility to radial highways
tax: full-value property-tax rate per USD 10,000
ptratio: pupil-teacher ratio by town
b: 1000(B - 0.63)^2 where B is the proportion of blacks by town
lstat: percentage of lower status of the population
medv: median value of owner-occupied homes in USD 1000’s
No lo introducimos en el código para no manchar mucho el Markdown, pero haciendo sapply(housing, unique), podemos comprobar qué valores toman cada una de las variables, y usamos esto para decidir qué variables son continuas, cuáles son discretas pero toman un rango suficientemente grande de valores y cuáles son categóricas. Vemos que la variable rad toma solo los valores 1,2,3,4,5,6,7,8,24. Siendo el rango de valores tan pequeño, nos vemos tentados de introducirla como variable categórica; probaremos de ambas formas, tanto de forma categórica como numérica. Chas es una variable categórica claramente, como viene especificado arriba. El resto de las variables son o bien continuas o discretas que toman suficientes valores como para dejarlas como numéricas. Vamos a ver, en caso de ser variable categórica, cuántas observaciones habría en cada grupo de la variable rad.
require(ggplot2)
ggplot(data = housing, aes(x = as.factor(rad))) +
geom_bar(fill = "red")
Vemos que la mayoría de las observaciones están en 4, 5 y 24. Sin embargo, inicialmente usaremos esta variable como numérica ya que en teoría es un índice de la accesibilidad a las “radial highways”; por lo tanto la diferencia numérica representa, a priori, una diferencia proporcional en la facilidad para acceder a estas carreteras, y por lo tanto podría tener una influencia en el precio mediano de las viviendas del distrito.
library(psych)##
## Attaching package: 'psych'## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha## The following object is masked from 'package:car':
##
## logitmulti.hist(housing[ , unlist(lapply(housing, is.numeric))])
library(fitdistrplus)
housing$chas <- as.factor(housing$chas)
for(i in 1:ncol(housing)) {
require(ggplot2)
if(is.numeric(housing[ , i])) {
hist(housing[ , i], main = names(housing)[i], col = "blue")
boxplot(housing[ , i], col = "green", main = names(housing)[i])
if(is.integer(housing[ , i])) {
descdist(housing[ , i], discrete = T) #si todos los elementos del vector son enteros, es una variable discreta que toma suficientes valores.
} else {
descdist(housing[ , i], discrete = F) #sino, es una variable continua
}
} else if(is.factor(housing[ , i])) {
counts <- table(housing[ , i])
barplot(counts, main = names(housing)[i])
descdist(as.numeric(housing[ , i]), discrete = T)
} else {
next
}
}
Para la variable crim, vemos que la mayoría de los valores se encuentran cerca del 0, con muuy pocas observaciones de distritos (towns) con un ratio de crimen per capita superior a 20. En el boxplot de esta misma variable nos queda claro que es una variable con una distribución asimétrica fuerte a la derecha, de hecho el square skewness es altísimo, como vemos en el tercer gráfico. La distribución que más se acerca a la de esta variable es una beta.
La variable zn tiene un comportamiento bastante similar, con muchas observaciones en valores cerca del 0 o 0, y pocas observaciones por encima de 20 “residential land zoned for lots” por cada 25,000 metros cuadrados. La distribución es asimétrica a la derecha, con una fuerte asimetría como puede apreciarse en el boxplot de esta variable. Esta variable tiene una distribución de tipo beta, según el gráfico de Cullen y Frey.
La variable indus nos muestra que la mayoría de “towns” (distritos) de Boston tienen una proporción de “non-retail business acres” entre 0 y 15, siendo esta porción de la distribución algo muy parecido a una Normal; sin embargo los valores más presentes están cerca del 20, con más de 150 observaciones; esto hace que la muestra pudiera dividirse casi en 2 partes en esta variable: aquellas observaciones con un indus de menos de 15 y las que tienen indus de más de 15. Como vemos en el boxplot, esta variable está más centrada y no tiene una asimetría tan clara como en los otros casos; ninguna observación, según el boxplot, sería considerada outlier en una normal. Curiosamente, vemos que la distribución que más se asemeja a esta variable es la uniforme.
Para la variable chas, podemos ver que la mayoría de los distritos no tienen cercanía al río Charles, ya que la categoría más abundante en esta variable categórica es, sin duda, el 0 (no cercanía al río Charles). Quizás nos faltan observaciones de la clase 1 para poder utilizar esta variable como predictor robusto, ya que las clases incluidas deben de ser suficientemente representativas.
La variable nox se distribuye de forma algo más parecida a la Normal; con un mínimo en 0.3850 y un máximo en 0.8710. Tiene una pequeña asimetría a la derecha, como podemos ver por la diferencia entre media y mediana. Sin embargo, esta asimetría no es tan fuerte como para las variables zn y crime, como vemos en el boxplot, ya que no quedan, a priori, variables que serían “outliers” en una Normal. La distribución que más se asemeja es también una beta, aunque con algo menos de asimetría, como vemos por el square of skewness del gráfico de Cullen y Frey. La curtosis también es algo inferior.
La variable rm tiene las observaciones más centradas que el resto de variables, y gráficamente vemos que se parece un poco más a la Normal que la mayoría del resto de variables que tenemos en el set de datos. En el boxplot vemos, sin embargo, que quedan bastantes observaciones fuera de los whiskers, lo cual nos indica que en una curva normal equivalente tendríamos varias observaciones como “outliers”; esto se puede deber a la gran centralización de los datos al rededor de la media. De hecho, como se puede ver en el gráfico de Cullen y Frey, se podría comportar como una Normal con una curtosis de 5 en lugar de 3.
En la variable age, vemos que la mayoría de los distritos tienen una proporción de viviendas ocupadas construidas antes de 1940 muy alta, en muchos caso casi del 100%. La distribución, si no fuera por este último rango de valores cerca del 100%, podría asemejarse bastante a una uniforme, ya que para el resto de valores de la variable las alturas de los “bins” son casi iguales (tenemos distritos con todos los tipos de proporciones de viviendas ocupadas construidas antes de 1940). Debido a este último de los “bins” observamos una asimetría a la izquierda en el boxplot. Como comentaba previamente, si no fuera por esta asimetría tendríamos una uniforme (gráfico de Cullen y Frey); pero la distribución que más se parece es una beta.
Para la variable dis, podemos ver que tiene una asimetría a la derecha a simple vista a juzgar por el histograma; esto lo confirmamos con el boxplot, que nos muestra que algunas observaciones quedarían como outliers en lo que sería una Normal. De nuevo vemos que la distribución que más se asemeja a la de esta variable sería una beta.
La variable rad, que podríamos haber dejado como categórica, vemos que tiene valores del 1 al 8 y luego el 24. Podríamos decir que, si no fuera por la diferencia de tamaños (para rad = 4, 5, 24 la altura o número de muestras es superior al resto), podría parecerse a una uniforme. De todas formas, es muy difícil atribuirle una distribución concreta a esta variable por contener números discretos de forma consecutiva hasta el 8 y luego el 24, que es 3 veces el anterior valor. Cualquier estimación de tendencia central o de cualquier otro momento del proceso estocástico que determina la distribución de rad estará sesgado por este hecho. Según el gráfico de Cullen y Frey esta variable podría seguir una distribución beta.
La variable tax tiene una distribuición muy parecida a rad, lo que nos hace pensar que muy probablemente estén relacionadas (puede ser que una de ellas determine a la otra). Las mismas consideraciones en cuanto a distribución posible etc que se hicieron para rad se podrían hacer para tax.
Para el ptratio observamos una asimetría a la izquierda, en el boxplot se ve que algunas observaciones quedarían como “outliers” en una Normal; la distribución más factible basándonos en el gráfico de Cullen y Frey sería una beta; tiene una curtosis algo inferior a la normal (3), con una asimetría mayor.
La variable b nos muestra que la mayoría de los distritos tienen un valor para la variable b cercano a 400. Casi todas las observaciones están en el rango más alto de la distribución, es decir tiene una clara asimetría a la izquierda. Para sacar el porcentaje exacto de negros por distrito habría que hacer: \(B = \frac{\sqrt{b}}{1000} + 0.63\) Vemos en el boxplot que hay muchas observaciones que serían outliers en una Normal, esto se debe al gran número de distritos con un b alto. La distribución que más se asemeja a b es una beta con una curtosis y asimetría cuadrada altísimos.
La variable lstat tiene una pequeña asimetría a la derecha, algo más suave que en otras variables. Algunas observaciones quedarían como “outliers” en una Normal, y la distribución que más se asemeja a lstat es una beta. Vemos, tanto en el histograma como en el boxplot, que la mayoría de los distritos de Boston tienen un porcentaje de gente de clase “baja” (“lower status”) bastante pequeño, de en torno al 10% (aunque algunos distritos superan el 30%).
Por último, la variable medv tiene también asimetría a la derecha; la mayoría de las casas tienen un valor mediano de entre 10,000 y 30,000$, aunque algunos distritos llegan a 50,000 de valor mediano de sus casas “ocupadas por propietarios”. En el boxplot vemos que algunos valores quedarían como “outliers” en una Normal. Por último, parece que la distribución que mejor se asemeja a cómo se distribuye medv es la beta, aunque también queda cerca la gamma.
Observando los gráficos de kurtosis-skewness, podemos ver que casi todas nuestras variables se distribuyen como una beta. Esto es positivo ya que, en general, cuando realizamos combinaciones lineales de variables que se distribuyen de la misma forma, estas combinaciones resultantes se distribuirán como una normal.
Antes de continuar con la exploración de variables, vamos a realizar tests de normalidad para la variable que queremos predecir, medv; además, dibujaremos un histograma y volveremos a realizar estos tests para 3 transformaciones diferentes de esta variable. En el Modelo Lineal se asume que las variable explicativas son determinísticas, y como \(Y = X\beta + \epsilon\), las constantes estimadas \(\hat{\beta}\) multiplicadas por la matriz “determinística” X, será determinística. De modo que, teóricamente, siguiendo esa ecuación, \(\epsilon\) e Y se distribuirán de igual forma. Como necesitamos que \(\epsilon \sim N(0, \sigma_\epsilon)\), estamos interesados en que la variable a predecir medv se comporte de forma lo más parecido posible a una Normal. Además, todos los test que usamos en el modelo para contrastar si los coeficientes estimados son significativos y demás, se basan en la suposición de normalidad. Por eso
shapiro.test(housing$medv)##
## Shapiro-Wilk normality test
##
## data: housing$medv
## W = 0.91718, p-value = 4.941e-16ad.test(housing$medv)##
## Anderson-Darling normality test
##
## data: housing$medv
## A = 11.822, p-value < 2.2e-16lillie.test(housing$medv)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: housing$medv
## D = 0.14919, p-value < 2.2e-16medv_sqrt <- sqrt(housing$medv)
hist(medv_sqrt, main = "historam for sqrt(medv)")
ad.test(medv_sqrt)##
## Anderson-Darling normality test
##
## data: medv_sqrt
## A = 4.9495, p-value = 2.954e-12lillie.test(medv_sqrt)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: medv_sqrt
## D = 0.11072, p-value = 2.703e-16shapiro.test(medv_sqrt)##
## Shapiro-Wilk normality test
##
## data: medv_sqrt
## W = 0.96988, p-value = 1.099e-08medv_log <- log(housing$medv)
hist(medv_log, main = "histogram for log(medv)")
shapiro.test(medv_log)##
## Shapiro-Wilk normality test
##
## data: medv_log
## W = 0.97574, p-value = 1.935e-07ad.test(medv_log)##
## Anderson-Darling normality test
##
## data: medv_log
## A = 3.83, p-value = 1.481e-09lillie.test(medv_log)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: medv_log
## D = 0.080923, p-value = 2.079e-08medv_logsq <- log(housing$medv)^2
hist(medv_logsq, main = "histogram for log squared medv")
shapiro.test(medv_logsq)##
## Shapiro-Wilk normality test
##
## data: medv_logsq
## W = 0.9787, p-value = 9.581e-07ad.test(medv_logsq)##
## Anderson-Darling normality test
##
## data: medv_logsq
## A = 3.7055, p-value = 2.964e-09lillie.test(medv_logsq)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: medv_logsq
## D = 0.098058, p-value = 1.305e-12Como vemos, algunas transformaciones posibles que probablemente darían mejor resultado en la regresión que la variable en bruto medv, serían hacer su raíz cuadrada o sacar el logaritmo natural.
require(GGally)
pairs(housing[ , unlist(lapply(housing, is.numeric))], col = "red")
ggcorr(housing, label = T)## Warning in ggcorr(housing, label = T): data in column(s) 'chas' are not
## numeric and were ignored
for (i in 1:13) {
plot(x = housing[ , i], y = housing[, "medv"], type = "p", col = "blue", main = paste0(names(housing)[i], " vs medv"), xlab = names(housing)[i], ylab = "medv")
}
Observando el gráfico de la matriz de correlaciones, vemos que vamos a tener problemas de multicolinearidad, ya que hay algunas variables que están muy correlacionadas; como nox e indus, tax y rad, nox y tax, indus y tax. Las correlaciones altas con la variable medv nos interesan, el problema es cuando \(X_j \approx \alpha_1*X_1 + \alpha_2*X_2 + ...+ \alpha_k*X_k\); es decir, cuando alguna de las variables independientes del modelo es aproximadamente una combinación lineal de otras variables independientes.
En los gráficos siguientes se ve que algunas variables tienen una clara relación con medv, mientras que en otros casos no está tan claro. Dejaremos primero que hablen los datos y los algoritmos de selección de variables, y posteriormente analizaremos una a una las variables introducidas, y describiremos los cambios o modificaciones que podemos hacerles para ganar poder predictivo.
2.B
Para comenzar a seleccionar variables del modelo antes de realizar transformaciones y demás, vamos a probar a hacer un step desde el modelo simple \(Y = \beta_0\) hasta el modelo en el que incluimos todas las variables excepto medv y todas las interacciones entre las mismas. Para hacer esto podemos basarnos en este link; al incluir las variables en la forma Y ~ X1 * X2 * X3 estamos expresando: \(Y = \beta_0 + \beta_1\cdot X1 + \beta_2 \cdot X2 + \beta_3 \cdot X3 + \beta_4 \cdot X1 \cdot X2 + \beta_5 \cdot X1 \cdot X3 + \beta_6 \cdot X2 \cdot X3 + \beta_7 \cdot X1 \cdot X2 \cdot X3 + \epsilon\). Como esto es computacionalmente muy costoso, y llevaría mucho tiempo encontrar la combinación óptima de variables en base al AIC, dejamos el fullmodel de este tipo comentado y utilizamos uno hecho a mano, con más restricciones, es decir sin comprobar todas las combinaciones posibles entre las variables.
library(leaps)
library(MASS)
nullmodel <- lm(medv ~ 1, data = housing)
#fullmodel <- lm(medv ~ crim*zn*indus*chas*nox*rm*age*dis*rad*tax*ptratio*b*lstat, data = housing)
fullmodel = lm(medv ~ crim + zn + indus + chas + nox + rm + age + dis + rad + tax + ptratio + b + lstat + crim:zn + crim:nox + crim:rm + crim*lstat + lstat:rad + lstat:b + lstat:rm + ptratio:nox + ptratio:lstat, data = housing)
step.both <- stepAIC(nullmodel, scope = list(lower = nullmodel, upper = fullmodel), direction = "both")## Start: AIC=2246.51
## medv ~ 1
##
## Df Sum of Sq RSS AIC
## + lstat 1 23243.9 19472 1851.0
## + rm 1 20654.4 22062 1914.2
## + ptratio 1 11014.3 31702 2097.6
## + indus 1 9995.2 32721 2113.6
## + tax 1 9377.3 33339 2123.1
## + nox 1 7800.1 34916 2146.5
## + crim 1 6440.8 36276 2165.8
## + rad 1 6221.1 36495 2168.9
## + age 1 6069.8 36647 2171.0
## + zn 1 5549.7 37167 2178.1
## + b 1 4749.9 37966 2188.9
## + dis 1 2668.2 40048 2215.9
## + chas 1 1312.1 41404 2232.7
## <none> 42716 2246.5
##
## Step: AIC=1851.01
## medv ~ lstat
##
## Df Sum of Sq RSS AIC
## + rm 1 4033.1 15439 1735.6
## + ptratio 1 2670.1 16802 1778.4
## + chas 1 786.3 18686 1832.2
## + dis 1 772.4 18700 1832.5
## + age 1 304.3 19168 1845.0
## + tax 1 274.4 19198 1845.8
## + b 1 198.3 19274 1847.8
## + zn 1 160.3 19312 1848.8
## + crim 1 146.9 19325 1849.2
## + indus 1 98.7 19374 1850.4
## <none> 19472 1851.0
## + rad 1 25.1 19447 1852.4
## + nox 1 4.8 19468 1852.9
## - lstat 1 23243.9 42716 2246.5
##
## Step: AIC=1735.58
## medv ~ lstat + rm
##
## Df Sum of Sq RSS AIC
## + rm:lstat 1 4343.6 11096 1570.4
## + ptratio 1 1711.3 13728 1678.1
## + chas 1 548.5 14891 1719.3
## + b 1 512.3 14927 1720.5
## + tax 1 425.2 15014 1723.5
## + dis 1 351.2 15088 1725.9
## + crim 1 311.4 15128 1727.3
## + rad 1 180.5 15259 1731.6
## + indus 1 61.1 15378 1735.6
## <none> 15439 1735.6
## + zn 1 56.6 15383 1735.7
## + age 1 20.2 15419 1736.9
## + nox 1 14.9 15424 1737.1
## - rm 1 4033.1 19472 1851.0
## - lstat 1 6622.6 22062 1914.2
##
## Step: AIC=1570.42
## medv ~ lstat + rm + lstat:rm
##
## Df Sum of Sq RSS AIC
## + ptratio 1 684.4 10411 1540.2
## + dis 1 466.8 10629 1550.7
## + crim 1 412.1 10684 1553.3
## + chas 1 365.3 10730 1555.5
## + age 1 172.3 10923 1564.5
## + b 1 106.3 10990 1567.5
## + tax 1 105.4 10990 1567.6
## <none> 11096 1570.4
## + zn 1 20.8 11075 1571.5
## + nox 1 16.5 11079 1571.7
## + rad 1 12.7 11083 1571.8
## + indus 1 7.7 11088 1572.1
## - lstat:rm 1 4343.6 15439 1735.6
##
## Step: AIC=1540.2
## medv ~ lstat + rm + ptratio + lstat:rm
##
## Df Sum of Sq RSS AIC
## + dis 1 561.2 9850.2 1514.2
## + chas 1 287.2 10124.2 1528.0
## + crim 1 252.9 10158.4 1529.8
## + age 1 215.3 10196.0 1531.6
## + zn 1 123.8 10287.5 1536.2
## + b 1 92.6 10318.8 1537.7
## + indus 1 57.5 10353.8 1539.4
## <none> 10411.4 1540.2
## + rad 1 38.6 10372.7 1540.3
## + nox 1 6.5 10404.9 1541.9
## + tax 1 3.7 10407.6 1542.0
## + ptratio:lstat 1 2.2 10409.1 1542.1
## - ptratio 1 684.4 11095.7 1570.4
## - lstat:rm 1 3316.6 13728.0 1678.1
##
## Step: AIC=1514.17
## medv ~ lstat + rm + ptratio + dis + lstat:rm
##
## Df Sum of Sq RSS AIC
## + crim 1 420.7 9429.5 1494.1
## + nox 1 349.8 9500.4 1497.9
## + chas 1 185.1 9665.1 1506.6
## + b 1 152.1 9698.1 1508.3
## + tax 1 119.7 9730.5 1510.0
## + indus 1 69.0 9781.2 1512.6
## <none> 9850.2 1514.2
## + zn 1 16.2 9834.0 1515.3
## + ptratio:lstat 1 4.9 9845.3 1515.9
## + rad 1 1.5 9848.7 1516.1
## + age 1 0.3 9849.9 1516.2
## - dis 1 561.2 10411.4 1540.2
## - ptratio 1 778.8 10629.0 1550.7
## - lstat:rm 1 3378.7 13228.9 1661.4
##
## Step: AIC=1494.08
## medv ~ lstat + rm + ptratio + dis + crim + lstat:rm
##
## Df Sum of Sq RSS AIC
## + nox 1 254.1 9175.3 1482.3
## + chas 1 154.6 9274.9 1487.7
## + crim:lstat 1 149.6 9279.8 1488.0
## + rad 1 96.1 9333.4 1490.9
## + b 1 54.4 9375.1 1493.2
## + zn 1 53.3 9376.2 1493.2
## + indus 1 44.5 9385.0 1493.7
## <none> 9429.5 1494.1
## + tax 1 9.1 9420.4 1495.6
## + ptratio:lstat 1 2.4 9427.0 1496.0
## + age 1 1.3 9428.2 1496.0
## + crim:rm 1 1.1 9428.3 1496.0
## - crim 1 420.7 9850.2 1514.2
## - ptratio 1 580.8 10010.2 1522.3
## - dis 1 729.0 10158.4 1529.8
## - lstat:rm 1 3565.9 12995.4 1654.4
##
## Step: AIC=1482.26
## medv ~ lstat + rm + ptratio + dis + crim + nox + lstat:rm
##
## Df Sum of Sq RSS AIC
## + rad 1 305.10 8870.2 1467.1
## + crim:lstat 1 195.81 8979.5 1473.3
## + chas 1 189.32 8986.0 1473.7
## + crim:nox 1 137.48 9037.9 1476.6
## + zn 1 55.83 9119.5 1481.2
## <none> 9175.3 1482.3
## + b 1 28.63 9146.7 1482.7
## + tax 1 19.07 9156.3 1483.2
## + age 1 9.61 9165.7 1483.7
## + nox:ptratio 1 9.31 9166.0 1483.7
## + ptratio:lstat 1 6.28 9169.1 1483.9
## + crim:rm 1 2.33 9173.0 1484.1
## + indus 1 0.01 9175.3 1484.3
## - nox 1 254.12 9429.5 1494.1
## - crim 1 325.04 9500.4 1497.9
## - ptratio 1 701.73 9877.1 1517.5
## - dis 1 943.41 10118.8 1529.8
## - lstat:rm 1 3152.56 12327.9 1629.7
##
## Step: AIC=1467.14
## medv ~ lstat + rm + ptratio + dis + crim + nox + rad + lstat:rm
##
## Df Sum of Sq RSS AIC
## + rad:lstat 1 725.2 8145.0 1426.0
## + tax 1 216.0 8654.3 1456.7
## + chas 1 182.2 8688.0 1458.6
## + crim:lstat 1 173.2 8697.0 1459.2
## + b 1 67.8 8802.5 1465.3
## + crim:nox 1 64.9 8805.3 1465.4
## + ptratio:lstat 1 49.0 8821.2 1466.3
## + nox:ptratio 1 43.0 8827.2 1466.7
## <none> 8870.2 1467.1
## + age 1 26.3 8844.0 1467.6
## + zn 1 26.3 8844.0 1467.6
## + crim:rm 1 22.6 8847.7 1467.9
## + indus 1 5.4 8864.8 1468.8
## - rad 1 305.1 9175.3 1482.3
## - nox 1 463.2 9333.4 1490.9
## - crim 1 580.1 9450.3 1497.2
## - ptratio 1 975.7 9845.9 1518.0
## - dis 1 1015.5 9885.8 1520.0
## - lstat:rm 1 3248.6 12118.8 1623.0
##
## Step: AIC=1425.98
## medv ~ lstat + rm + ptratio + dis + crim + nox + rad + lstat:rm +
## lstat:rad
##
## Df Sum of Sq RSS AIC
## + tax 1 256.2 7888.8 1411.8
## + chas 1 83.5 8061.5 1422.8
## + zn 1 67.0 8078.0 1423.8
## + indus 1 61.4 8083.6 1424.2
## <none> 8145.0 1426.0
## + crim:lstat 1 20.4 8124.6 1426.7
## + b 1 20.1 8124.9 1426.7
## + ptratio:lstat 1 17.2 8127.8 1426.9
## + crim:rm 1 14.5 8130.5 1427.1
## + nox:ptratio 1 9.4 8135.7 1427.4
## + crim:nox 1 0.2 8144.8 1428.0
## + age 1 0.0 8145.0 1428.0
## - crim 1 281.5 8426.5 1441.2
## - nox 1 645.7 8790.7 1462.6
## - lstat:rad 1 725.2 8870.2 1467.1
## - dis 1 906.1 9051.1 1477.4
## - ptratio 1 1086.0 9231.1 1487.3
## - lstat:rm 1 3876.8 12021.8 1621.0
##
## Step: AIC=1411.81
## medv ~ lstat + rm + ptratio + dis + crim + nox + rad + tax +
## lstat:rm + lstat:rad
##
## Df Sum of Sq RSS AIC
## + zn 1 138.3 7750.5 1404.9
## + ptratio:lstat 1 67.4 7821.4 1409.5
## + chas 1 58.2 7830.6 1410.1
## <none> 7888.8 1411.8
## + b 1 13.3 7875.5 1413.0
## + crim:lstat 1 11.0 7877.8 1413.1
## + crim:rm 1 7.7 7881.1 1413.3
## + crim:nox 1 4.1 7884.7 1413.5
## + indus 1 1.5 7887.3 1413.7
## + age 1 0.0 7888.8 1413.8
## + nox:ptratio 1 0.0 7888.8 1413.8
## - tax 1 256.2 8145.0 1426.0
## - crim 1 275.6 8164.4 1427.2
## - nox 1 463.2 8352.0 1438.7
## - lstat:rad 1 765.5 8654.3 1456.7
## - dis 1 897.5 8786.3 1464.3
## - ptratio 1 987.1 8875.9 1469.5
## - lstat:rm 1 3864.2 11753.0 1611.5
##
## Step: AIC=1404.86
## medv ~ lstat + rm + ptratio + dis + crim + nox + rad + tax +
## zn + lstat:rm + lstat:rad
##
## Df Sum of Sq RSS AIC
## + chas 1 56.4 7694.1 1403.2
## + ptratio:lstat 1 40.0 7710.5 1404.2
## <none> 7750.5 1404.9
## + crim:zn 1 23.3 7727.2 1405.3
## + crim:lstat 1 15.9 7734.6 1405.8
## + b 1 13.3 7737.2 1406.0
## + crim:nox 1 7.4 7743.1 1406.4
## + crim:rm 1 7.4 7743.1 1406.4
## + nox:ptratio 1 4.1 7746.4 1406.6
## + age 1 0.8 7749.7 1406.8
## + indus 1 0.1 7750.4 1406.9
## - zn 1 138.3 7888.8 1411.8
## - crim 1 292.0 8042.5 1421.6
## - tax 1 327.5 8078.0 1423.8
## - nox 1 425.8 8176.3 1429.9
## - ptratio 1 691.7 8442.2 1446.1
## - lstat:rad 1 835.9 8586.4 1454.7
## - dis 1 1009.9 8760.4 1464.8
## - lstat:rm 1 3676.5 11427.0 1599.3
##
## Step: AIC=1403.17
## medv ~ lstat + rm + ptratio + dis + crim + nox + rad + tax +
## zn + chas + lstat:rm + lstat:rad
##
## Df Sum of Sq RSS AIC
## + ptratio:lstat 1 40.6 7653.5 1402.5
## + crim:zn 1 34.2 7659.9 1402.9
## <none> 7694.1 1403.2
## + crim:lstat 1 13.7 7680.4 1404.3
## + b 1 12.4 7681.7 1404.3
## + crim:nox 1 10.8 7683.3 1404.5
## + crim:rm 1 8.1 7686.0 1404.6
## + nox:ptratio 1 4.9 7689.2 1404.8
## - chas 1 56.4 7750.5 1404.9
## + indus 1 0.8 7693.3 1405.1
## + age 1 0.4 7693.7 1405.1
## - zn 1 136.5 7830.6 1410.1
## - crim 1 285.3 7979.5 1419.6
## - tax 1 298.9 7993.0 1420.5
## - nox 1 443.6 8137.8 1429.5
## - ptratio 1 661.7 8355.8 1442.9
## - lstat:rad 1 739.5 8433.6 1447.6
## - dis 1 983.1 8677.2 1462.0
## - lstat:rm 1 3533.6 11227.7 1592.4
##
## Step: AIC=1402.49
## medv ~ lstat + rm + ptratio + dis + crim + nox + rad + tax +
## zn + chas + lstat:rm + lstat:rad + lstat:ptratio
##
## Df Sum of Sq RSS AIC
## <none> 7653.5 1402.5
## + nox:ptratio 1 28.8 7624.7 1402.6
## + crim:zn 1 26.6 7626.9 1402.7
## - lstat:ptratio 1 40.6 7694.1 1403.2
## + crim:lstat 1 13.5 7640.0 1403.6
## + b 1 13.4 7640.1 1403.6
## + crim:nox 1 11.6 7641.9 1403.7
## + crim:rm 1 7.2 7646.3 1404.0
## - chas 1 57.0 7710.5 1404.2
## + indus 1 1.8 7651.7 1404.4
## + age 1 0.0 7653.5 1404.5
## - zn 1 109.1 7762.6 1407.7
## - crim 1 280.3 7933.8 1418.7
## - nox 1 316.6 7970.1 1421.0
## - tax 1 332.0 7985.5 1422.0
## - lstat:rad 1 754.8 8408.3 1448.1
## - dis 1 792.7 8446.2 1450.3
## - lstat:rm 1 3422.3 11075.8 1587.5summary(step.both)##
## Call:
## lm(formula = medv ~ lstat + rm + ptratio + dis + crim + nox +
## rad + tax + zn + chas + lstat:rm + lstat:rad + lstat:ptratio,
## data = housing)
##
## Residuals:
## Min 1Q Median 3Q Max
## -25.4893 -2.0974 -0.3293 1.5893 25.2042
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.069776 5.945622 0.853 0.39425
## lstat 2.048506 0.370269 5.532 5.14e-08 ***
## rm 8.799980 0.494855 17.783 < 2e-16 ***
## ptratio -1.011371 0.215731 -4.688 3.58e-06 ***
## dis -1.164229 0.163096 -7.138 3.41e-12 ***
## crim -0.119037 0.028043 -4.245 2.62e-05 ***
## nox -14.224407 3.152994 -4.511 8.06e-06 ***
## rad 0.744310 0.083356 8.929 < 2e-16 ***
## tax -0.013300 0.002879 -4.620 4.90e-06 ***
## zn 0.030797 0.011629 2.648 0.00835 **
## chas1 1.383941 0.722703 1.915 0.05608 .
## lstat:rm -0.479545 0.032331 -14.832 < 2e-16 ***
## lstat:rad -0.026843 0.003854 -6.966 1.05e-11 ***
## lstat:ptratio 0.026059 0.016128 1.616 0.10679
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.944 on 492 degrees of freedom
## Multiple R-squared: 0.8208, Adjusted R-squared: 0.8161
## F-statistic: 173.4 on 13 and 492 DF, p-value: < 2.2e-16Ahora vamos a ver si encontramos un codo para restringir el número de variables con este último modelo que nos da el step.
model <- regsubsets( medv ~ lstat + rm + ptratio + dis + crim + nox +
rad + tax + zn + chas + lstat:rm + lstat:rad + lstat:ptratio, data = housing, nvmax = 12)
plot(model)
plot(summary(model)$rsq, type = "l")
No se termina de ver un codo claro, aunque hay pequeños cambios de pendiente del R^2 en 2, 3 y 8 variables, por lo que nos quedaremos de momento con lstat, rm, ptratio, dis, crim, nox, rad (tax la eliminamos porque está muy correlacionada con rad, contienen información redundante), zn, lstat:rm y lstat:rad.
require(car)
crPlots(lm(medv ~ lstat + rm + ptratio + dis + crim + nox +
rad + zn , data = housing))
Basándonos en los gráficos que vemos arriba, realizamos algunas transformaciones en determinadas variables (ver código) para “linearizar” las relaciones entre las mismas y así mejorar el modelo.
housing2 <- housing
housing2$dissq <- housing2$dis^2
housing2$medv_sqrt <- sqrt(housing2$medv)
modelo_1 <- lm( medv_sqrt ~ log(lstat) + rm + I(rm^2) + log(ptratio) + log(dissq) + sqrt(crim) + nox +
log(rad) + log(ptratio) + lstat:rad + lstat:I(rm^2), data = housing2)
summary(modelo_1)##
## Call:
## lm(formula = medv_sqrt ~ log(lstat) + rm + I(rm^2) + log(ptratio) +
## log(dissq) + sqrt(crim) + nox + log(rad) + log(ptratio) +
## lstat:rad + lstat:I(rm^2), data = housing2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.14393 -0.22451 -0.01328 0.21531 2.16690
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.3536214 1.0982513 13.980 < 2e-16 ***
## log(lstat) -0.4952239 0.1384393 -3.577 0.000381 ***
## rm -1.6084125 0.2744394 -5.861 8.42e-09 ***
## I(rm^2) 0.1534358 0.0204996 7.485 3.30e-13 ***
## log(ptratio) -1.2770306 0.1703709 -7.496 3.07e-13 ***
## log(dissq) -0.2572464 0.0311758 -8.251 1.43e-15 ***
## sqrt(crim) -0.1576843 0.0263640 -5.981 4.25e-09 ***
## nox -2.0321305 0.3066601 -6.627 8.99e-11 ***
## log(rad) 0.2754266 0.0403328 6.829 2.52e-11 ***
## lstat:rad -0.0008881 0.0002548 -3.486 0.000535 ***
## I(rm^2):lstat -0.0006155 0.0003344 -1.840 0.066306 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3915 on 495 degrees of freedom
## Multiple R-squared: 0.8286, Adjusted R-squared: 0.8251
## F-statistic: 239.3 on 10 and 495 DF, p-value: < 2.2e-16plot(modelo_1)
vif(modelo_1)## log(lstat) rm I(rm^2) log(ptratio) log(dissq)
## 22.799173 122.500646 114.137191 1.489964 3.728735
## sqrt(crim) nox log(rad) lstat:rad I(rm^2):lstat
## 4.823076 4.160277 4.101816 8.263450 17.596038En primer lugar, viendo los p-value, tenemos que tomar una decisión de qué hacemos con \(rm^2 \cdot lstat\), ya que si usamos un \(\alpha = 0.05\), tendríamos que eliminarla. La cuestión es la siguiente: en caso de que la variable mencionada sea realmente relevante (significativa), si la eliminamos pasaremos esta variable al error, incumpliendo supuestos y por lo tanto sesgando el estimador (ya no tendríamos un estimador lineal insesgado según el teorema de Gauss-Markov); si resulta no ser relevante y la incluimos en el modelo, dejaremos de tener un BLUE (Best Linear Unbiased Estimator), pero seguiremos teniendo un estimador insesgado. Ante este tradeoff, tomo la decisión, debido a la cercanía del p-value con la frontera (\(\alpha\)) de rechazo, y a esta explicación teórica previa, de mantener la variable en el modelo. Es preferible perder precisión que sesgar el estimador.
Este modelo tiene un R^2 ajustao bastante alto, de 0.8251, pero la colinearidad, basados en el Variance Inflation Factor, está presente en muchas de las variables que hemos introducido. A partir de este modelo iremos tratando de afinar un poco más, eliminando algunas variables que puedan resultar redundantes para lidiar con el problema de la colinearidad.
En cuanto a posibles observaciones influyentes (outliers que sesgan el estimador), sólo encontramos una observación en la primera frontera de la distancia de Cook de 0.5; probamos a eliminarla y ver qué tal se ajusta el modelo en ese caso.
modelo_1_no_highleverage <- lm( medv_sqrt ~ log(lstat) + rm + I(rm^2) + log(ptratio) + log(dissq) + sqrt(crim) + nox +
log(rad) + log(ptratio) + lstat:rad + lstat:I(rm^2), data = housing2[-365, ])
summary(modelo_1_no_highleverage)##
## Call:
## lm(formula = medv_sqrt ~ log(lstat) + rm + I(rm^2) + log(ptratio) +
## log(dissq) + sqrt(crim) + nox + log(rad) + log(ptratio) +
## lstat:rad + lstat:I(rm^2), data = housing2[-365, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.3702 -0.2237 -0.0162 0.2044 2.1729
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.6723572 1.0616708 14.762 < 2e-16 ***
## log(lstat) -0.4287606 0.1341110 -3.197 0.00148 **
## rm -1.9747376 0.2717339 -7.267 1.44e-12 ***
## I(rm^2) 0.1863332 0.0205184 9.081 < 2e-16 ***
## log(ptratio) -1.1333997 0.1661819 -6.820 2.66e-11 ***
## log(dissq) -0.2429051 0.0301929 -8.045 6.45e-15 ***
## sqrt(crim) -0.1573408 0.0254549 -6.181 1.33e-09 ***
## nox -1.8615952 0.2974092 -6.259 8.40e-10 ***
## log(rad) 0.2922364 0.0390399 7.486 3.29e-13 ***
## lstat:rad -0.0010236 0.0002470 -4.144 4.02e-05 ***
## I(rm^2):lstat -0.0007063 0.0003232 -2.185 0.02934 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.378 on 494 degrees of freedom
## Multiple R-squared: 0.8405, Adjusted R-squared: 0.8373
## F-statistic: 260.4 on 10 and 494 DF, p-value: < 2.2e-16Como podemos ver, la eliminación de esta observación ha resultado en un considerable aumento de R^2 ajustado, además de un aumento en el valor absoluto del t-value de la variable \(lstat \cdot rm^2\), lo cual se traduce en que resulta significativa para \(\alpha = 0.5\).
modelo_2_2 <- lm(log(medv) ~ log(crim) + log(nox) + age + rm + dis + rad + tax + ptratio + log(b) + log(lstat), data = housing)
summary(modelo_2_2)##
## Call:
## lm(formula = log(medv) ~ log(crim) + log(nox) + age + rm + dis +
## rad + tax + ptratio + log(b) + log(lstat), data = housing)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.96367 -0.09945 -0.00030 0.10485 0.87582
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.8538906 0.2150498 17.921 < 2e-16 ***
## log(crim) -0.0251261 0.0112267 -2.238 0.025660 *
## log(nox) -0.2922704 0.1035332 -2.823 0.004950 **
## age 0.0009901 0.0005782 1.713 0.087422 .
## rm 0.0665823 0.0183769 3.623 0.000321 ***
## dis -0.0398687 0.0077819 -5.123 4.31e-07 ***
## rad 0.0109240 0.0030343 3.600 0.000350 ***
## tax -0.0005858 0.0001391 -4.210 3.03e-05 ***
## ptratio -0.0314861 0.0051006 -6.173 1.40e-09 ***
## log(b) 0.0589428 0.0127915 4.608 5.18e-06 ***
## log(lstat) -0.4113530 0.0267254 -15.392 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2007 on 495 degrees of freedom
## Multiple R-squared: 0.7637, Adjusted R-squared: 0.7589
## F-statistic: 160 on 10 and 495 DF, p-value: < 2.2e-16Este modelo empeora significativamente con respecto al anterior, en términos de R^2 ajustado.
Este chunk de aquí abajo se irá utilizando para probar todos los modelos que se nos ocurran; sin embargo, para no alargar este Markdown más de la cuenta, se presentarán sólo los modelos definitivos, con variaciones para poder compararles en el Ejercicio 3.
library(car)
modelo_2_3 <- lm(sqrt(medv) ~ sqrt(crim) + nox + rm + I(rm^2) + log(dis) + log(rad) + log(ptratio) + sqrt(b) + log(lstat), data = housing2)
summary(modelo_2_3)##
## Call:
## lm(formula = sqrt(medv) ~ sqrt(crim) + nox + rm + I(rm^2) + log(dis) +
## log(rad) + log(ptratio) + sqrt(b) + log(lstat), data = housing2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.94083 -0.22462 -0.00882 0.20747 2.26759
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.793732 0.947431 16.670 < 2e-16 ***
## sqrt(crim) -0.202898 0.022901 -8.860 < 2e-16 ***
## nox -1.887521 0.309275 -6.103 2.10e-09 ***
## rm -1.572882 0.256882 -6.123 1.87e-09 ***
## I(rm^2) 0.143527 0.020118 7.134 3.47e-12 ***
## log(dis) -0.510179 0.063029 -8.094 4.48e-15 ***
## log(rad) 0.204383 0.034705 5.889 7.17e-09 ***
## log(ptratio) -1.328173 0.172026 -7.721 6.43e-14 ***
## sqrt(b) 0.018844 0.005664 3.327 0.000943 ***
## log(lstat) -0.826996 0.049674 -16.648 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.396 on 496 degrees of freedom
## Multiple R-squared: 0.8243, Adjusted R-squared: 0.8211
## F-statistic: 258.6 on 9 and 496 DF, p-value: < 2.2e-16plot(modelo_2_3)
plot(hatvalues(modelo_2_3))
vif(modelo_2_3)## sqrt(crim) nox rm I(rm^2) log(dis)
## 3.557944 4.136872 104.926703 107.466428 3.725012
## log(rad) log(ptratio) sqrt(b) log(lstat)
## 2.969078 1.485076 1.337814 2.869665crPlots(modelo_2_3)
En principio, vemos que con las variables incluidas podemos explicar, teóricamente, el 82.11% de la varianza de \(\sqrt{medv}\).
Podríamos quitarle algunos puntos influyentes como la observación 365; sin embargo, dado que ningún outlier supera la distancia de Cook de 0.5, creo que es una buena decisión no eliminar estas variables, ya que estaríamos en típico caso de overfitting. Nos interesa que el modelo generalice bien a la hora de predecir, y si ajustamos perfectamente los coeficientes para conseguir un \(R^2\) alto en la fase de train, estaremos lastrando la capacidad del modelo de no quedar demasiado lejos de observaciones algo más extremas que la mayoría; esto nos hace perder algo de precisión en la muestra de datos que le estamos enseñando (en este caso es la muestra entera; comprobaremos la fiabilidad del modelo viendo diferentes datos en el siguiente ejercicio), pero es preferible a hacer overfitting en cuanto a procurarle robustez al modelo. Sin embargo, debemos mencionar que las observaciones 365, 372 y 373 son claramente outliers en el dataset mostrado según el Q-Q Plot. Testamos la normalidad de los residuos:
shapiro.test(modelo_2_3$residuals)##
## Shapiro-Wilk normality test
##
## data: modelo_2_3$residuals
## W = 0.94072, p-value = 2.528e-13shapiro.test(modelo_2_3$residuals[-c(365, 372, 373)])##
## Shapiro-Wilk normality test
##
## data: modelo_2_3$residuals[-c(365, 372, 373)]
## W = 0.98461, p-value = 3.647e-05shapiro.test(modelo_2_3$residuals[-c(365, 369, 372, 373)])##
## Shapiro-Wilk normality test
##
## data: modelo_2_3$residuals[-c(365, 369, 372, 373)]
## W = 0.98722, p-value = 0.0002201shapiro.test(modelo_2_3$residuals[-which(hatvalues(modelo_2_3) > 0.15)])##
## Shapiro-Wilk normality test
##
## data: modelo_2_3$residuals[-which(hatvalues(modelo_2_3) > 0.15)]
## W = 0.94103, p-value = 2.982e-13hist(modelo_2_3$residuals, main = "residuals histogram", xlab = "residuals", col = "blue")
Como podemos comprobar, \(p-val < \alpha, \alpha = 0.05\), de modo que rechazamos la hipótesis de que los residuos se comportan de forma normal. Esto puede suponer un problema para prácticamente cualquier estimador lineal, en caso de estar suficientemente lejos de la normal (lo cual por el bajísimo p-value es bastante probable). Todos los tests que realizamos los hacemos sobre la suposición de que \(\epsilon \sim N(0, \sigma_\epsilon)\), en cuanto esto no se cumple no podemos estar seguros de tener un buen modelo predictivo.
Vemos que al eliminar tan sólo esas 3 observaciones que tenían unos errores anormalmente grandes (siempre en relación a la muestra de datos) el p-value sube hasta \(3.647 \cdot 10^-5\). Y cuando quitamos además la observación 369, sube hasta \(0.0002201\). En función de si estas observaciones fueran verdaderos outliers o simplemente son una clase de casas que no está bien representada en esta muestra o para la que faltan variables explicativas tomaríamos una decisión u otra respecto a su eliminación. Para comprobar esto necesitaríamos más datos.
Vemos con el Variance Inflation Factor (VIF) que podríamos tener colinealidad entre \(rm\) y \(rm^2\), esto era posible ya que, pese a no ser esta segunda una combinación lineal perfecta de la primera, es la misma variable elevada al cuadrado, y por lo tanto existe la posibilidad de que contenga prácticamente la misma información que la primera.
Mantendremos además del modelo modelo_2_3, uno nuevo en el que no incluiremos rm, tan sólo su transformación cuadrática \(rm^2\):
housing2$rmsq <- housing2$rm^2
modelo_no_colineal <- lm(sqrt(medv) ~ sqrt(crim) + nox + rm + log(dissq) + log(rad) + log(ptratio) + sqrt(b) + log(lstat), data = housing2)
summary(modelo_no_colineal)##
## Call:
## lm(formula = sqrt(medv) ~ sqrt(crim) + nox + rm + log(dissq) +
## log(rad) + log(ptratio) + sqrt(b) + log(lstat), data = housing2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.26123 -0.24472 -0.02506 0.22205 2.11067
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.008327 0.701893 15.684 < 2e-16 ***
## sqrt(crim) -0.184704 0.023874 -7.737 5.74e-14 ***
## nox -2.036284 0.323691 -6.291 6.93e-10 ***
## rm 0.242938 0.036453 6.664 7.07e-11 ***
## log(dissq) -0.302458 0.032317 -9.359 < 2e-16 ***
## log(rad) 0.211355 0.036391 5.808 1.13e-08 ***
## log(ptratio) -1.540591 0.177732 -8.668 < 2e-16 ***
## sqrt(b) 0.021392 0.005929 3.608 0.00034 ***
## log(lstat) -0.893091 0.051194 -17.445 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4154 on 497 degrees of freedom
## Multiple R-squared: 0.8063, Adjusted R-squared: 0.8032
## F-statistic: 258.6 on 8 and 497 DF, p-value: < 2.2e-16crPlots(modelo_no_colineal)
otro_modelo_alternativo <- lm(log(medv) ~ sqrt(crim) + nox + rmsq + log(dis) + log(rad) + log(ptratio) + sqrt(b) + log(lstat), data = housing2)
summary(otro_modelo_alternativo)##
## Call:
## lm(formula = log(medv) ~ sqrt(crim) + nox + rmsq + log(dis) +
## log(rad) + log(ptratio) + sqrt(b) + log(lstat), data = housing2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.80359 -0.09703 -0.00378 0.10177 0.80611
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.772722 0.290188 19.893 < 2e-16 ***
## sqrt(crim) -0.111612 0.010651 -10.479 < 2e-16 ***
## nox -0.756601 0.144542 -5.234 2.45e-07 ***
## rmsq 0.006850 0.001274 5.376 1.18e-07 ***
## log(dis) -0.218492 0.028933 -7.552 2.08e-13 ***
## log(rad) 0.101482 0.016250 6.245 9.10e-10 ***
## log(ptratio) -0.582304 0.079646 -7.311 1.07e-12 ***
## sqrt(b) 0.009339 0.002644 3.532 0.000451 ***
## log(lstat) -0.364858 0.023126 -15.777 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1854 on 497 degrees of freedom
## Multiple R-squared: 0.7975, Adjusted R-squared: 0.7942
## F-statistic: 244.7 on 8 and 497 DF, p-value: < 2.2e-16Con este último modelo alternativo configuramos los cuatro modelos que vamos a usar para el ejercicio 3:
Modelo_2_3 (\(R_adj^2\) = 0.8211)
Modelo_no_colineal (\(R_adj^2\) = 0.8032)
Otro_modelo_alternativo (\(R_adj^2\) = 0.7918)
Modelo_1 (\(R_adj^2\) = 0.8251)
La homocedasticidad y normalidad de los residuos del mejor modelo serán estudiados también en el siguiente ejercicio. De todas las modificaciones que le hemos hecho al primer modelo que sacamos con stepAIC, las que mejor han resultado en términos de \(R^2_adj\) son las introducidas en el Modelo_1. Por lo tanto, de momento seleccionamos ese modelo como el definitivo. Utilizaremos los resultados del próximo ejercicio para salir de dudas.
EJERCICIO 3.
Todos los comentarios realizados en el anterior ejercicio al respecto de los supuestos, de las decisiones de incluir o no variables, de tratar de no sesgar al estimador y hacerle ganar precisión etc, serán comprobados en este ejercicio. En definitiva, el cross validation nos sirve para realizar un proceso iterativo y robusto de train y test splits para comprobar la capacidad predictiva real de los modelos, entrenados y estimados sus parámetros con una porción de los datos y comprobando qué tal predice ese modelo unos nuevos datos que no ha visto.
Esto se hará el número de observaciones de la muestra para el Leave One Out Cross Validation y 10 veces para el 10-Fold Cross Validation. Después, para elegir el modelo definitivo, nos podríamos guiar por distintos criterios. En concreto el output de las funciones de caret nos brinda el RMSE, R^2 y MAE de la predicción. Si la variable dependiente en todos los modelos estuviera exactamente en las mismas unidades, sería una buena idea utilizar el RMSE; sin embargo, al ser en algunos casos \(\sqrt{medv}\) y en otros \(ln(medv)\), utilizaremos como medida para elegir el modelo definitivo el R^2, que nos da el acierto en la predicción en forma de porcentaje, de 0 a 1, de tal forma que podemos comparar el performance de los modelos.
require(caret)## Loading required package: caret## Loading required package: lattice##
## Attaching package: 'caret'## The following object is masked from 'package:survival':
##
## clusterLEAVE ONE OUT CROSS VALIDATION.
set.seed(1)
train_control <- trainControl(method="LOOCV")
model1 <- train(sqrt(medv) ~ sqrt(crim) + nox + rm + I(rm^2) + log(dis) + log(rad) + log(ptratio) + sqrt(b) + log(lstat), data= housing2, trControl=train_control, method="lm")
rmse1 <- as.numeric(model1$results["RMSE"])
rsq1 <- as.numeric(model1$results["Rsquared"])
model1_best <- model1$finalModel
summary(model1)##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.94083 -0.22462 -0.00882 0.20747 2.26759
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.793732 0.947431 16.670 < 2e-16 ***
## `sqrt(crim)` -0.202898 0.022901 -8.860 < 2e-16 ***
## nox -1.887521 0.309275 -6.103 2.10e-09 ***
## rm -1.572882 0.256882 -6.123 1.87e-09 ***
## `I(rm^2)` 0.143527 0.020118 7.134 3.47e-12 ***
## `log(dis)` -0.510179 0.063029 -8.094 4.48e-15 ***
## `log(rad)` 0.204383 0.034705 5.889 7.17e-09 ***
## `log(ptratio)` -1.328173 0.172026 -7.721 6.43e-14 ***
## `sqrt(b)` 0.018844 0.005664 3.327 0.000943 ***
## `log(lstat)` -0.826996 0.049674 -16.648 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.396 on 496 degrees of freedom
## Multiple R-squared: 0.8243, Adjusted R-squared: 0.8211
## F-statistic: 258.6 on 9 and 496 DF, p-value: < 2.2e-16print(model1)## Linear Regression
##
## 506 samples
## 8 predictor
##
## No pre-processing
## Resampling: Leave-One-Out Cross-Validation
## Summary of sample sizes: 505, 505, 505, 505, 505, 505, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 0.4041713 0.8133138 0.2880238
##
## Tuning parameter 'intercept' was held constant at a value of TRUEmodel2 <- train(sqrt(medv) ~ sqrt(crim) + nox + rm + log(dissq) + log(rad) + log(ptratio) + sqrt(b) + log(lstat), data = housing2, trControl = train_control, method = "lm")
summary(model2)##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.26123 -0.24472 -0.02506 0.22205 2.11067
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.008327 0.701893 15.684 < 2e-16 ***
## `sqrt(crim)` -0.184704 0.023874 -7.737 5.74e-14 ***
## nox -2.036284 0.323691 -6.291 6.93e-10 ***
## rm 0.242938 0.036453 6.664 7.07e-11 ***
## `log(dissq)` -0.302458 0.032317 -9.359 < 2e-16 ***
## `log(rad)` 0.211355 0.036391 5.808 1.13e-08 ***
## `log(ptratio)` -1.540591 0.177732 -8.668 < 2e-16 ***
## `sqrt(b)` 0.021392 0.005929 3.608 0.00034 ***
## `log(lstat)` -0.893091 0.051194 -17.445 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4154 on 497 degrees of freedom
## Multiple R-squared: 0.8063, Adjusted R-squared: 0.8032
## F-statistic: 258.6 on 8 and 497 DF, p-value: < 2.2e-16print(model2)## Linear Regression
##
## 506 samples
## 8 predictor
##
## No pre-processing
## Resampling: Leave-One-Out Cross-Validation
## Summary of sample sizes: 505, 505, 505, 505, 505, 505, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 0.4218599 0.7965889 0.3126714
##
## Tuning parameter 'intercept' was held constant at a value of TRUErmse2 <- as.numeric(model2$results["RMSE"])
rsq2 <- as.numeric(model2$results["Rsquared"])
model2_best <- model2$finalModel
model3 <- train(log(medv) ~ sqrt(crim) + nox + I(rm^2) + log(dis) + log(rad) + log(ptratio) + sqrt(b) + log(lstat), data = housing2, trControl = train_control, method = "lm")
summary(model3)##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.80359 -0.09703 -0.00378 0.10177 0.80611
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.772722 0.290188 19.893 < 2e-16 ***
## `sqrt(crim)` -0.111612 0.010651 -10.479 < 2e-16 ***
## nox -0.756601 0.144542 -5.234 2.45e-07 ***
## `I(rm^2)` 0.006850 0.001274 5.376 1.18e-07 ***
## `log(dis)` -0.218492 0.028933 -7.552 2.08e-13 ***
## `log(rad)` 0.101482 0.016250 6.245 9.10e-10 ***
## `log(ptratio)` -0.582304 0.079646 -7.311 1.07e-12 ***
## `sqrt(b)` 0.009339 0.002644 3.532 0.000451 ***
## `log(lstat)` -0.364858 0.023126 -15.777 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1854 on 497 degrees of freedom
## Multiple R-squared: 0.7975, Adjusted R-squared: 0.7942
## F-statistic: 244.7 on 8 and 497 DF, p-value: < 2.2e-16print(model3)## Linear Regression
##
## 506 samples
## 8 predictor
##
## No pre-processing
## Resampling: Leave-One-Out Cross-Validation
## Summary of sample sizes: 505, 505, 505, 505, 505, 505, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 0.188469 0.7870218 0.1348858
##
## Tuning parameter 'intercept' was held constant at a value of TRUErmse3 <- as.numeric(model3$results["RMSE"])
rsq3 <- as.numeric(model3$results["Rsquared"])
model3_best <- model3$finalModel
model4 <- train(medv_sqrt ~ log(lstat) + rm + I(rm^2) + log(ptratio) + log(dissq) + sqrt(crim) + nox +
log(rad) + log(ptratio) + lstat:rad + lstat:I(rm^2), data = housing2, trControl = train_control, method = "lm")
summary(model4)##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.14393 -0.22451 -0.01328 0.21531 2.16690
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.3536214 1.0982513 13.980 < 2e-16 ***
## `log(lstat)` -0.4952239 0.1384393 -3.577 0.000381 ***
## rm -1.6084125 0.2744394 -5.861 8.42e-09 ***
## `I(rm^2)` 0.1534358 0.0204996 7.485 3.30e-13 ***
## `log(ptratio)` -1.2770306 0.1703709 -7.496 3.07e-13 ***
## `log(dissq)` -0.2572464 0.0311758 -8.251 1.43e-15 ***
## `sqrt(crim)` -0.1576843 0.0263640 -5.981 4.25e-09 ***
## nox -2.0321305 0.3066601 -6.627 8.99e-11 ***
## `log(rad)` 0.2754266 0.0403328 6.829 2.52e-11 ***
## `lstat:rad` -0.0008881 0.0002548 -3.486 0.000535 ***
## `I(rm^2):lstat` -0.0006155 0.0003344 -1.840 0.066306 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3915 on 495 degrees of freedom
## Multiple R-squared: 0.8286, Adjusted R-squared: 0.8251
## F-statistic: 239.3 on 10 and 495 DF, p-value: < 2.2e-16print(model4)## Linear Regression
##
## 506 samples
## 7 predictor
##
## No pre-processing
## Resampling: Leave-One-Out Cross-Validation
## Summary of sample sizes: 505, 505, 505, 505, 505, 505, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 0.4022555 0.8151037 0.2874553
##
## Tuning parameter 'intercept' was held constant at a value of TRUErmse4 <- as.numeric(model4$results["RMSE"])
rsq4 <- as.numeric(model4$results["Rsquared"])
model4_best <- model4$finalModel
rmse_loocv <- as.vector(c(rmse1, rmse2, rmse3, rmse4))
rsqvec <- as.vector(c(rsq1, rsq2, rsq3, rsq4))
#print(paste("el mejor rmse es ", min(rmse_loocv),"que pertenece al modelo", which(rmse_loocv == min(rmse_loocv))) )
print(paste("el mejor r^2 es ", max(rsqvec), "que pertenece al modelo", which(rsqvec == max(rsqvec))))## [1] "el mejor r^2 es 0.815103696805514 que pertenece al modelo 4"Por el método de Leave One Out CV sería una buena idea elegir el modelo número 4, es decir, Modelo_1. Tiene, de media, un R^2 de 0.8151 en los test set (la media del R^2 en cada test set de los 506 sets). Testamos la homocedasticidad y normalidad de los residuos de este modelo.
require(lmtest)
shapiro.test(model4_best$residuals)##
## Shapiro-Wilk normality test
##
## data: model4_best$residuals
## W = 0.93942, p-value = 1.716e-13bptest(model4_best)##
## studentized Breusch-Pagan test
##
## data: model4_best
## BP = 56.203, df = 10, p-value = 1.881e-08Sin embargo, por buenas que sean las predicciones, seguimos rechazando tanto la hipótesis de normalidad como la de homocedasticidad. Por eso debemos tener cuidado, y siempre tener en mente que por buenas que sean las medidas de error que estemos obteniendo es posible que los test que realizamos y la confianza en el modelo no estén justificadas por incumplimiento de los supuestos del modelo.
Comprobaremos si el 10-Fold Cross Validation, que en teoría debería ser más robusto, nos indica que elijamos ese modelo o por el contrario nos obliga a tomar una decisión difícil.
10-FOLD CROSS VALIDATION
set.seed(1)
train_control <- trainControl(method="cv", number=10)
model1 <- train(sqrt(medv) ~ sqrt(crim) + nox + rm + I(rm^2) + log(dis) + log(rad) + log(ptratio) + sqrt(b) + log(lstat), data= housing2, trControl=train_control, method="lm")
rmse1 <- as.numeric(model1$results["RMSE"])
rsq1 <- as.numeric(model1$results["Rsquared"])
model1_best <- model1$finalModel
summary(model1)##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.94083 -0.22462 -0.00882 0.20747 2.26759
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.793732 0.947431 16.670 < 2e-16 ***
## `sqrt(crim)` -0.202898 0.022901 -8.860 < 2e-16 ***
## nox -1.887521 0.309275 -6.103 2.10e-09 ***
## rm -1.572882 0.256882 -6.123 1.87e-09 ***
## `I(rm^2)` 0.143527 0.020118 7.134 3.47e-12 ***
## `log(dis)` -0.510179 0.063029 -8.094 4.48e-15 ***
## `log(rad)` 0.204383 0.034705 5.889 7.17e-09 ***
## `log(ptratio)` -1.328173 0.172026 -7.721 6.43e-14 ***
## `sqrt(b)` 0.018844 0.005664 3.327 0.000943 ***
## `log(lstat)` -0.826996 0.049674 -16.648 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.396 on 496 degrees of freedom
## Multiple R-squared: 0.8243, Adjusted R-squared: 0.8211
## F-statistic: 258.6 on 9 and 496 DF, p-value: < 2.2e-16print(model2)## Linear Regression
##
## 506 samples
## 8 predictor
##
## No pre-processing
## Resampling: Leave-One-Out Cross-Validation
## Summary of sample sizes: 505, 505, 505, 505, 505, 505, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 0.4218599 0.7965889 0.3126714
##
## Tuning parameter 'intercept' was held constant at a value of TRUEmodel2 <- train(sqrt(medv) ~ sqrt(crim) + nox + rm + log(dissq) + log(rad) + log(ptratio) + sqrt(b) + log(lstat), data = housing2, trControl = train_control, method = "lm")
summary(model2)##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.26123 -0.24472 -0.02506 0.22205 2.11067
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.008327 0.701893 15.684 < 2e-16 ***
## `sqrt(crim)` -0.184704 0.023874 -7.737 5.74e-14 ***
## nox -2.036284 0.323691 -6.291 6.93e-10 ***
## rm 0.242938 0.036453 6.664 7.07e-11 ***
## `log(dissq)` -0.302458 0.032317 -9.359 < 2e-16 ***
## `log(rad)` 0.211355 0.036391 5.808 1.13e-08 ***
## `log(ptratio)` -1.540591 0.177732 -8.668 < 2e-16 ***
## `sqrt(b)` 0.021392 0.005929 3.608 0.00034 ***
## `log(lstat)` -0.893091 0.051194 -17.445 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4154 on 497 degrees of freedom
## Multiple R-squared: 0.8063, Adjusted R-squared: 0.8032
## F-statistic: 258.6 on 8 and 497 DF, p-value: < 2.2e-16print(model2)## Linear Regression
##
## 506 samples
## 8 predictor
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 455, 455, 457, 456, 456, 454, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 0.4218683 0.8031355 0.3161113
##
## Tuning parameter 'intercept' was held constant at a value of TRUErmse2 <- as.numeric(model2$results["RMSE"])
rsq2 <- as.numeric(model2$results["Rsquared"])
model2_best <- model2$finalModel
model3 <- train(log(medv) ~ sqrt(crim) + nox + I(rm^2) + log(dis) + log(rad) + log(ptratio) + sqrt(b) + log(lstat), data = housing2, trControl = train_control, method = "lm")
summary(model3)##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.80359 -0.09703 -0.00378 0.10177 0.80611
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.772722 0.290188 19.893 < 2e-16 ***
## `sqrt(crim)` -0.111612 0.010651 -10.479 < 2e-16 ***
## nox -0.756601 0.144542 -5.234 2.45e-07 ***
## `I(rm^2)` 0.006850 0.001274 5.376 1.18e-07 ***
## `log(dis)` -0.218492 0.028933 -7.552 2.08e-13 ***
## `log(rad)` 0.101482 0.016250 6.245 9.10e-10 ***
## `log(ptratio)` -0.582304 0.079646 -7.311 1.07e-12 ***
## `sqrt(b)` 0.009339 0.002644 3.532 0.000451 ***
## `log(lstat)` -0.364858 0.023126 -15.777 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1854 on 497 degrees of freedom
## Multiple R-squared: 0.7975, Adjusted R-squared: 0.7942
## F-statistic: 244.7 on 8 and 497 DF, p-value: < 2.2e-16print(model3)## Linear Regression
##
## 506 samples
## 8 predictor
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 455, 455, 456, 457, 455, 456, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 0.1890429 0.7889447 0.1358331
##
## Tuning parameter 'intercept' was held constant at a value of TRUErmse3 <- as.numeric(model3$results["RMSE"])
rsq3 <- as.numeric(model3$results["Rsquared"])
model3_best <- model3$finalModel
model4 <- train(medv_sqrt ~ log(lstat) + rm + I(rm^2) + log(ptratio) + log(dissq) + sqrt(crim) + nox +
log(rad) + log(ptratio) + lstat:rad + lstat:I(rm^2), data = housing2, trControl = train_control, method = "lm")
summary(model4)##
## Call:
## lm(formula = .outcome ~ ., data = dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.14393 -0.22451 -0.01328 0.21531 2.16690
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.3536214 1.0982513 13.980 < 2e-16 ***
## `log(lstat)` -0.4952239 0.1384393 -3.577 0.000381 ***
## rm -1.6084125 0.2744394 -5.861 8.42e-09 ***
## `I(rm^2)` 0.1534358 0.0204996 7.485 3.30e-13 ***
## `log(ptratio)` -1.2770306 0.1703709 -7.496 3.07e-13 ***
## `log(dissq)` -0.2572464 0.0311758 -8.251 1.43e-15 ***
## `sqrt(crim)` -0.1576843 0.0263640 -5.981 4.25e-09 ***
## nox -2.0321305 0.3066601 -6.627 8.99e-11 ***
## `log(rad)` 0.2754266 0.0403328 6.829 2.52e-11 ***
## `lstat:rad` -0.0008881 0.0002548 -3.486 0.000535 ***
## `I(rm^2):lstat` -0.0006155 0.0003344 -1.840 0.066306 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3915 on 495 degrees of freedom
## Multiple R-squared: 0.8286, Adjusted R-squared: 0.8251
## F-statistic: 239.3 on 10 and 495 DF, p-value: < 2.2e-16print(model4)## Linear Regression
##
## 506 samples
## 7 predictor
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 455, 455, 456, 456, 456, 456, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 0.3988127 0.8233312 0.288492
##
## Tuning parameter 'intercept' was held constant at a value of TRUErmse4 <- as.numeric(model4$results["RMSE"])
rsq4 <- as.numeric(model4$results["Rsquared"])
model4_best <- model4$finalModel
rmse_10fcv <- as.vector(c(rmse1, rmse2, rmse3, rmse4))
rsqvec <- as.vector(c(rsq1, rsq2, rsq3, rsq4))
print(paste("el mejor r^2 es ", max(rsqvec), "que pertenece al modelo", which(rsqvec == max(rsqvec))))## [1] "el mejor r^2 es 0.823331150029626 que pertenece al modelo 4"Este método nos indica que el modelo que mejor predice de media en las 10 muestras aleatorias diferentes presentadas tras estimar sus \(\beta\), es el modelo número 4; de tal forma que este modelo es el que mejor predice basándonos en estos dos métodos de cross validation.
Por lo tanto, determinamos que el modelo que mejor predice es el modelo número 4, es decir Modelo_1:medv_sqrt ~ log(lstat) + rm + I(rm^2) + log(ptratio) + log(dissq) + sqrt(crim) + nox + log(rad) + log(ptratio) + lstat:rad + lstat:I(rm^2). La decisión provisional del ejercicio 2 de quedarnos con este modelo como el definitivo resulta ser validada por ambos métodos de validación cruzada, ya que además de ser el que mejor bondad de ajuste tenía viendo todos los datos, es el que mejor predice de media en múltiples particiones train-test.
Como la homocedasticidad y la normalidad ya fueron rechazadas para este mismo modelo en el LOOCV, concluiremos diciendo que decidimos quedarnos definitivamente con este modelo pues tras hacer cross validation ha sido el que mejor predice de media; sin embargo es un modelo aún con muchas carencias, ya que no se cumplen algunos de los supuestos del Teorema de Gauss-Markov y eso provoca que nuestro estimador pueda estar sesgado, además de no ser el más eficiente. El problema de la heterocedasticidad podría resolverse haciendo la transformación de las variables, ajustándolas por la varianza del error, de tal forma que el modelo resultante tenga unos errores homocedásticos, siguiendo este ejemplo.